home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 7 / Apprentice-Release7.iso / Environments / PowerLisp 2.01 / Supplemental Documentation / Documentation / Chapter 12. Numbers < prev    next >
Text File  |  1995-03-27  |  132KB  |  3,162 lines
  1. Common Lisp the Language, 2nd Edition
  2. -------------------------------------------------------------------------------
  3.  
  4. 12. Numbers
  5.  
  6. Common Lisp provides several different representations for numbers. These
  7. representations may be divided into four categories: integers, ratios,
  8. floating-point numbers, and complex numbers. Many numeric functions will accept
  9. any kind of number; they are generic. Other functions accept only certain kinds
  10. of numbers.
  11.  
  12. [change_begin]
  13. Note that this remark, predating the design of the Common Lisp Object System,
  14. uses the term ``generic'' in a generic sense and not necessarily in the
  15. technical sense used by CLOS (see chapter 2).
  16. [change_end]
  17.  
  18. In general, numbers in Common Lisp are not true objects; eq cannot be counted
  19. upon to operate on them reliably. In particular, it is possible that the
  20. expression
  21.  
  22. (let ((x z) (y z)) (eq x y))
  23.  
  24. may be false rather than true if the value of z is a number.
  25.  
  26. -------------------------------------------------------------------------------
  27. Rationale: This odd breakdown of eq in the case of numbers allows the
  28. implementor enough design freedom to produce exceptionally efficient numerical
  29. code on conventional architectures. MacLisp requires this freedom, for example,
  30. in order to produce compiled numerical code equal in speed to Fortran. Common
  31. Lisp makes this same restriction, if not for this freedom, then at least for
  32. the sake of compatibility.
  33. -------------------------------------------------------------------------------
  34.  
  35. If two objects are to be compared for ``identity,'' but either might be a
  36. number, then the predicate eql is probably appropriate; if both objects are
  37. known to be numbers, then = may be preferable.
  38.  
  39. -------------------------------------------------------------------------------
  40.  
  41.    *  Precision, Contagion, and Coercion
  42.    *  Predicates on Numbers
  43.    *  Comparisons on Numbers
  44.    *  Arithmetic Operations
  45.    *  Irrational and Transcendental Functions
  46.         o  Exponential and Logarithmic Functions
  47.         o  Trigonometric and Related Functions
  48.         o  Branch Cuts, Principal Values, and Boundary Conditions in the
  49.           Complex Plane
  50.    *  Type Conversions and Component Extractions on Numbers
  51.    *  Logical Operations on Numbers
  52.    *  Byte Manipulation Functions
  53.    *  Random Numbers
  54.    *  Implementation Parameters
  55.  
  56. -------------------------------------------------------------------------------
  57.  
  58. 12.1. Precision, Contagion, and Coercion
  59.  
  60. In general, computations with floating-point numbers are only approximate. The
  61. precision of a floating-point number is not necessarily correlated at all with
  62. the accuracy of that number. For instance, 3.142857142857142857 is a more
  63. precise approximation to than 3.14159, but the latter is more accurate. The
  64. precision refers to the number of bits retained in the representation. When an
  65. operation combines a short floating-point number with a long one, the result
  66. will be a long floating-point number. This rule is made to ensure that as much
  67. accuracy as possible is preserved; however, it is by no means a guarantee.
  68. Common Lisp numerical routines do assume, however, that the accuracy of an
  69. argument does not exceed its precision. Therefore when two small floating-point
  70. numbers are combined, the result will always be a small floating-point number.
  71. This assumption can be overridden by first explicitly converting a small
  72. floating-point number to a larger representation. (Common Lisp never converts
  73. automatically from a larger size to a smaller one.)
  74.  
  75. Rational computations cannot overflow in the usual sense (though of course
  76. there may not be enough storage to represent one), as integers and ratios may
  77. in principle be of any magnitude. Floating-point computations may get exponent
  78. overflow or underflow; this is an error.
  79.  
  80. [change_begin]
  81. X3J13 voted in June 1989 (FLOAT-UNDERFLOW)   to address certain problems
  82. relating to floating-point overflow and underflow, but certain parts of the
  83. proposed solution were not adopted, namely to add the macro
  84. without-floating-underflow-traps to the language and to require certain
  85. behavior of floating-point overflow and underflow. The committee agreed that
  86. this area of the language requires more discussion before a solution is
  87. standardized.
  88.  
  89. For the record, the proposal that was considered and rejected (for the nonce)
  90. introduced a macro without-floating-underflow-traps that would execute its body
  91. in such a way that, within its dynamic extent, a floating-point underflow must
  92. not signal an error but instead must produce either a denormalized number or
  93. zero as the result. The rejected proposal also specified the following
  94. treatment of overflow and underflow:
  95.  
  96.    *  A floating-point computation that overflows should signal an error of
  97.      type floating-point-overflow.
  98.  
  99.    *  Unless the dynamic extent of a use of without-floating-underflow-traps, a
  100.      floating-point computation that underflows should signal an error of type
  101.      floating-point-underflow. A result that can be represented only in
  102.      denormalized form must be considered an underflow in implementations that
  103.      support denormalized floating-point numbers.
  104.  
  105. These points refer to conditions floating-point-overflow and
  106. floating-point-underflow that were approved by X3J13 and are described in
  107. section 29.5.
  108. [change_end]
  109.  
  110. When rational and floating-point numbers are compared or combined by a
  111. numerical function, the rule of floating-point contagion is followed: when a
  112. rational meets a floating-point number, the rational is first converted to a
  113. floating-point number of the same format. For functions such as + that take
  114. more than two arguments, it may be that part of the operation is carried out
  115. exactly using rationals and then the rest is done using floating-point
  116. arithmetic.
  117.  
  118. [change_begin]
  119. X3J13 voted in January 1989 (CONTAGION-ON-NUMERICAL-COMPARISONS)   to apply the
  120. rule of floating-point contagion stated above to the case of combining rational
  121. and floating-point numbers. For comparing, the following rule is to be used
  122. instead: When a rational number and a floating-point number are to be compared
  123. by a numerical function, in effect the floating-point number is first converted
  124. to a rational number as if by the function rational, and then an exact
  125. comparison of two rational numbers is performed. It is of course valid to use a
  126. more efficient implementation than actually calling the function rational, as
  127. long as the result of the comparison is the same. In the case of complex
  128. numbers, the real and imaginary parts are handled separately.
  129.  
  130. -------------------------------------------------------------------------------
  131. Rationale: In general, accuracy cannot be preserved in combining operations,
  132. but it can be preserved in comparisons, and preserving it makes that part of
  133. Common Lisp algebraically a bit more tractable. In particular, this change
  134. prevents the breakdown of transitivity. Let a be the result of (/ 10.0
  135. single-float-epsilon), and let j be the result of (floor a). (Note that (= a (+
  136. a 1.0)) is true, by the definition of single-float-epsilon.) Under the old
  137. rules, all of (<= a j), (< j (+ j 1)), and (<= (+ j 1) a) would be true;
  138. transitivity would then imply that (< a a) ought to be true, but of course it
  139. is false, and therefore transitivity fails. Under the new rule, however, (<= (+
  140. j 1) a) is false.
  141. -------------------------------------------------------------------------------
  142.  
  143. [change_end]
  144.  
  145. For functions that are mathematically associative (and possibly commutative), a
  146. Common Lisp implementation may process the arguments in any manner consistent
  147. with associative (and possibly commutative) rearrangement. This does not affect
  148. the order in which the argument forms are evaluated, of course; that order is
  149. always left to right, as in all Common Lisp function calls. What is left loose
  150. is the order in which the argument values are processed. The point of all this
  151. is that implementations may differ in which automatic coercions are applied
  152. because of differing orders of argument processing. As an example, consider
  153. this expression:
  154.  
  155. (+ 1/3 2/3 1.0D0 1.0 1.0E-15)
  156.  
  157. One implementation might process the arguments from left to right, first adding
  158. 1/3 and 2/3 to get 1, then converting that to a double-precision floating-point
  159. number for combination with 1.0D0, then successively converting and adding 1.0
  160. and 1.0E-15. Another implementation might process the arguments from right to
  161. left, first performing a single-precision floating-point addition of 1.0 and
  162. 1.0E-15 (and probably losing some accuracy in the process!), then converting
  163. the sum to double precision and adding 1.0D0, then converting 2/3 to
  164. double-precision floating-point and adding it, and then converting 1/3 and
  165. adding that. A third implementation might first scan all the arguments, process
  166. all the rationals first to keep that part of the computation exact, then find
  167. an argument of the largest floating-point format among all the arguments and
  168. add that, and then add in all other arguments, converting each in turn (all in
  169. a perhaps misguided attempt to make the computation as accurate as possible).
  170. In any case, all three strategies are legitimate. The user can of course
  171. control the order of processing explicitly by writing several calls; for
  172. example:
  173.  
  174. (+ (+ 1/3 2/3) (+ 1.0D0 1.0E-15) 1.0)
  175.  
  176. The user can also control all coercions simply by writing calls to coercion
  177. functions explicitly.
  178.  
  179. In general, then, the type of the result of a numerical function is a
  180. floating-point number of the largest format among all the floating-point
  181. arguments to the function; but if the arguments are all rational, then the
  182. result is rational (except for functions that can produce mathematically
  183. irrational results, in which case a single-format floating-point number may
  184. result).
  185.  
  186. There is a separate rule of complex contagion. As a rule, complex numbers never
  187. result from a numerical function unless one or more of the arguments is
  188. complex. (Exceptions to this rule occur among the irrational and transcendental
  189. functions, specifically expt, log, sqrt, asin, acos, acosh, and atanh; see
  190. section 12.5.) When a non-complex number meets a complex number, the
  191. non-complex number is in effect first converted to a complex number by
  192. providing an imaginary part of zero.
  193.  
  194. If any computation produces a result that is a ratio of two integers such that
  195. the denominator evenly divides the numerator, then the result is immediately
  196. converted to the equivalent integer. This is called the rule of rational
  197. canonicalization.
  198.  
  199. If the result of any computation would be a complex rational with a zero
  200. imaginary part, the result is immediately converted to a non-complex rational
  201. number by taking the real part. This is called the rule of complex
  202. canonicalization. Note that this rule does not apply to complex numbers whose
  203. components are floating-point numbers. Whereas #C(5 0) and 5 are not distinct
  204. values in Common Lisp (they are always eql), #C(5.0 0.0) and 5.0 are always
  205. distinct values in Common Lisp (they are never eql, although they are equalp).
  206.  
  207. -------------------------------------------------------------------------------
  208.  
  209. 12.2. Predicates on Numbers
  210.  
  211. Each of the following functions tests a single number for a specific property.
  212. Each function requires that its argument be a number; to call one with a
  213. non-number is an error.
  214.  
  215. [Function]
  216. zerop number
  217.  
  218. This predicate is true if number is zero (the integer zero, a floating-point
  219. zero, or a complex zero), and is false otherwise. Regardless of whether an
  220. implementation provides distinct representations for positive and negative
  221. floating-point zeros, (zerop -0.0) is always true. It is an error if the
  222. argument number is not a number.
  223.  
  224. [Function]
  225. plusp number
  226.  
  227. This predicate is true if number is strictly greater than zero, and is false
  228. otherwise. It is an error if the argument number is not a non-complex number.
  229.  
  230. [Function]
  231. minusp number
  232.  
  233. This predicate is true if number is strictly less than zero, and is false
  234. otherwise. Regardless of whether an implementation provides distinct
  235. representations for positive and negative floating-point zeros, (minusp -0.0)
  236. is always false. (The function float-sign may be used to distinguish a negative
  237. zero.) It is an error if the argument number is not a non-complex number.
  238.  
  239. [Function]
  240. oddp integer
  241.  
  242. This predicate is true if the argument integer is odd (not divisible by 2), and
  243. otherwise is false. It is an error if the argument is not an integer.
  244.  
  245. [Function]
  246. evenp integer
  247.  
  248. This predicate is true if the argument integer is even (divisible by 2), and
  249. otherwise is false. It is an error if the argument is not an integer.
  250.  
  251. See also the data-type predicates integerp, rationalp, floatp, complexp, and
  252. numberp.
  253.  
  254. -------------------------------------------------------------------------------
  255.  
  256. 12.3. Comparisons on Numbers
  257.  
  258. Each of the functions in this section requires that its arguments all be
  259. numbers; to call one with a non-number is an error. Unless otherwise specified,
  260. each works on all types of numbers, automatically performing any required
  261. coercions when arguments are of different types.
  262.  
  263. [Function]
  264. = number &rest more-numbers
  265. /= number &rest more-numbers
  266. < number &rest more-numbers
  267. > number &rest more-numbers
  268. <= number &rest more-numbers
  269. >= number &rest more-numbers
  270.  
  271. These functions each take one or more arguments. If the sequence of arguments
  272. satisfies a certain condition:
  273.  
  274. =            all the same
  275. /=           all different
  276. <            monotonically increasing
  277. >            monotonically decreasing
  278. <=           monotonically nondecreasing
  279. >=           monotonically nonincreasing
  280.  
  281. then the predicate is true, and otherwise is false. Complex numbers may be
  282. compared using = and /=, but the others require non-complex arguments. Two
  283. complex numbers are considered equal by = if their real parts are equal and
  284. their imaginary parts are equal according to =. A complex number may be
  285. compared with a non-complex number with = or /=. For example:
  286.  
  287. (= 3 3) is true.                (/= 3 3) is false.
  288. (= 3 5) is false.               (/= 3 5) is true.
  289. (= 3 3 3 3) is true.            (/= 3 3 3 3) is false.
  290. (= 3 3 5 3) is false.           (/= 3 3 5 3) is false.
  291. (= 3 6 5 2) is false.           (/= 3 6 5 2) is true.
  292. (= 3 2 3) is false.             (/= 3 2 3) is false.
  293. (< 3 5) is true.                (<= 3 5) is true.
  294. (< 3 -5) is false.              (<= 3 -5) is false.
  295. (< 3 3) is false.               (<= 3 3) is true.
  296. (< 0 3 4 6 7) is true.          (<= 0 3 4 6 7) is true.
  297. (< 0 3 4 4 6) is false.         (<= 0 3 4 4 6) is true.
  298. (> 4 3) is true.                (>= 4 3) is true.
  299. (> 4 3 2 1 0) is true.          (>= 4 3 2 1 0) is true.
  300. (> 4 3 3 2 0) is false.         (>= 4 3 3 2 0) is true.
  301. (> 4 3 1 2 0) is false.         (>= 4 3 1 2 0) is false.
  302. (= 3) is true.                  (/= 3) is true.
  303. (< 3) is true.                  (<= 3) is true.
  304. (= 3.0 #C(3.0 0.0)) is true.    (/= 3.0 #C(3.0 1.0)) is true.
  305. (= 3 3.0) is true.              (= 3.0s0 3.0d0) is true.
  306. (= 0.0 -0.0) is true.           (= 5/2 2.5) is true.
  307. (> 0.0 -0.0) is false.          (= 0 -0.0) is true.
  308.  
  309. With two arguments, these functions perform the usual arithmetic comparison
  310. tests. With three or more arguments, they are useful for range checks, as shown
  311. in the following example:
  312.  
  313. (<= 0 x 9)                      ;true if x is between 0 and 9, inclusive
  314. (< 0.0 x 1.0)                   ;true if x is between 0.0 and 1.0, exclusive
  315. (< -1 j (length s))             ;true if j is a valid index for s
  316. (<= 0 j k (- (length s) 1))     ;true if j and k are each valid
  317.                                 ; indices for s and jk
  318.  
  319. -------------------------------------------------------------------------------
  320. Rationale: The ``unequality'' relation is called /= rather than <> (the name
  321. used in Pascal) for two reasons. First, /= of more than two arguments is not
  322. the same as the or of < and > of those same arguments. Second, unequality is
  323. meaningful for complex numbers even though < and > are not. For both reasons it
  324. would be misleading to associate unequality with the names of < and >.
  325. -------------------------------------------------------------------------------
  326. Compatibility note: In Common Lisp, the comparison operations perform
  327. ``mixed-mode'' comparisons: (= 3 3.0) is true. In MacLisp, there must be
  328. exactly two arguments, and they must be either both fixnums or both
  329. floating-point numbers. To compare two numbers for numerical equality and type
  330. equality, use eql.
  331. -------------------------------------------------------------------------------
  332.  
  333. [Function]
  334. max number &rest more-numbers
  335. min number &rest more-numbers
  336.  
  337. The arguments may be any non-complex numbers. max returns the argument that is
  338. greatest (closest to positive infinity). min returns the argument that is least
  339. (closest to negative infinity).
  340.  
  341. For max, if the arguments are a mixture of rationals and floating-point
  342. numbers, and the largest argument is a rational, then the implementation is
  343. free to produce either that rational or its floating-point approximation; if
  344. the largest argument is a floating-point number of a smaller format than the
  345. largest format of any floating-point argument, then the implementation is free
  346. to return the argument in its given format or expanded to the larger format.
  347. More concisely, the implementation has the choice of returning the largest
  348. argument as is or applying the rules of floating-point contagion, taking all
  349. the arguments into consideration for contagion purposes. Also, if two or more
  350. of the arguments are equal, then any one of them may be chosen as the value to
  351. return. Similar remarks apply to min (replacing ``largest argument'' by
  352. ``smallest argument'').
  353.  
  354. (max 6 12) => 12                (min 6 12) => 6
  355. (max -6 -12) => -6              (min -6 -12) => -12
  356. (max 1 3 2 -7) => 3             (min 1 3 2 -7) => -7
  357. (max -2 3 0 7) => 7             (min -2 3 0 7) => -2
  358. (max 3) => 3                    (min 3) => 3
  359. (max 5.0 2) => 5.0              (min 5.0 2) => 2 or 2.0
  360. (max 3.0 7 1) => 7 or 7.0       (min 3.0 7 1) => 1 or 1.0
  361. (max 1.0s0 7.0d0) => 7.0d0
  362. (min 1.0s0 7.0d0) => 1.0s0 or 1.0d0
  363. (max 3 1 1.0s0 1.0d0) => 3 or 3.0d0
  364. (min 3 1 1.0s0 1.0d0) => 1 or 1.0s0 or 1.0d0
  365.  
  366. -------------------------------------------------------------------------------
  367.  
  368. 12.4. Arithmetic Operations
  369.  
  370. Each of the functions in this section requires that its arguments all be
  371. numbers; to call one with a non-number is an error. Unless otherwise specified,
  372. each works on all types of numbers, automatically performing any required
  373. coercions when arguments are of different types.
  374.  
  375. [Function]
  376. + &rest numbers
  377.  
  378. This returns the sum of the arguments. If there are no arguments, the result is
  379. 0, which is an identity for this operation.
  380.  
  381. -------------------------------------------------------------------------------
  382. Compatibility note: While + is compatible with its use in Lisp Machine Lisp, it
  383. is incompatible with MacLisp, which uses + for fixnum-only addition.
  384. -------------------------------------------------------------------------------
  385.  
  386. [Function]
  387. - number &rest more-numbers
  388.  
  389. The function -, when given one argument, returns the negative of that argument.
  390.  
  391. The function -, when given more than one argument, successively subtracts from
  392. the first argument all the others, and returns the result. For example, (- 3 4
  393. 5) => -6.
  394.  
  395. -------------------------------------------------------------------------------
  396. Compatibility note: While - is compatible with its use in Lisp Machine Lisp, it
  397. is incompatible with MacLisp, which uses - for fixnum-only subtraction. Also, -
  398. differs from difference as used in most Lisp systems in the case of one
  399. argument.
  400. -------------------------------------------------------------------------------
  401.  
  402. [Function]
  403. * &rest numbers
  404.  
  405. This returns the product of the arguments. If there are no arguments, the
  406. result is 1, which is an identity for this operation.
  407.  
  408. -------------------------------------------------------------------------------
  409. Compatibility note: While * is compatible with its use in Lisp Machine Lisp, it
  410. is incompatible with MacLisp, which uses * for fixnum-only multiplication.
  411. -------------------------------------------------------------------------------
  412.  
  413. [Function]
  414. / number &rest more-numbers
  415.  
  416. The function /, when given more than one argument, successively divides the
  417. first argument by all the others and returns the result.
  418.  
  419. [change_begin]
  420. It is generally accepted that it is an error for any argument other than the
  421. first to be zero.
  422. [change_end]
  423.  
  424. With one argument, / reciprocates the argument.
  425.  
  426. [change_begin]
  427. It is generally accepted that it is an error in this case for the argument to
  428. be zero.
  429. [change_end]
  430.  
  431. / will produce a ratio if the mathematical quotient of two integers is not an
  432. exact integer. For example:
  433.  
  434. (/ 12 4) => 3
  435. (/ 13 4) => 13/4
  436. (/ -8) => -1/8
  437. (/ 3 4 5) => 3/20
  438.  
  439. To divide one integer by another producing an integer result, use one of the
  440. functions floor, ceiling, truncate, or round.
  441.  
  442. If any argument is a floating-point number, then the rules of floating-point
  443. contagion apply.
  444.  
  445. -------------------------------------------------------------------------------
  446. Compatibility note: What / does is totally unlike what the usual // or quotient
  447. operator does. In most Lisp systems, quotient behaves like / except when
  448. dividing integers, in which case it behaves like truncate of two arguments;
  449. this behavior is mathematically intractable, leading to such anomalies as
  450.  
  451. (quotient 1.0 2.0) => 0.5   but   (quotient 1 2) => 0
  452.  
  453. In contrast, the Common Lisp function / produces these results:
  454.  
  455. (/ 1.0 2.0) => 0.5          and   (/ 1 2) => 1/2
  456.  
  457. In practice quotient is used only when one is sure that both arguments are
  458. integers, or when one is sure that at least one argument is a floating-point
  459. number. / is tractable for its purpose and works for any numbers.
  460. -------------------------------------------------------------------------------
  461.  
  462. [Function]
  463. 1+ number
  464. 1- number
  465.  
  466. (1+ x) is the same as (+ x 1).
  467.  
  468. (1- x) is the same as (- x 1). Note that the short name may be confusing: (1-
  469. x) does not mean 1-x; rather, it means x-1.
  470.  
  471. -------------------------------------------------------------------------------
  472. Rationale: These are included primarily for compatibility with MacLisp and Lisp
  473. Machine Lisp. Some programmers prefer always to write (+ x 1) and (- x 1)
  474. instead of (1+ x) and (1- x).
  475. -------------------------------------------------------------------------------
  476. Implementation note: Compiler writers are very strongly encouraged to ensure
  477. that (1+ x) and (+ x 1) compile into identical code, and similarly for (1- x)
  478. and (- x 1), to avoid pressure on a Lisp programmer to write possibly less
  479. clear code for the sake of efficiency. This can easily be done as a
  480. source-language transformation.
  481. -------------------------------------------------------------------------------
  482.  
  483. [Macro]
  484. incf place [delta]
  485. decf place [delta]
  486.  
  487. The number produced by the form delta is added to (incf) or subtracted from
  488. (decf) the number in the generalized variable named by place, and the sum is
  489. stored back into place and returned. The form place may be any form acceptable
  490. as a generalized variable to setf. If delta is not supplied, then the number in
  491. place is changed by 1. For example:
  492.  
  493. (setq n 0)
  494. (incf n) => 1                   and now n => 1
  495. (decf n 3) => -2                and now n => -2
  496. (decf n -5) => 3                and now n => 3
  497. (decf n) => 2                   and now n => 2
  498.  
  499. The effect of (incf place delta) is roughly equivalent to
  500.  
  501. (setf place (+ place delta))
  502.  
  503. except that the latter would evaluate any subforms of place twice, whereas incf
  504. takes care to evaluate them only once. Moreover, for certain place forms incf
  505. may be significantly more efficient than the setf version.
  506.  
  507. [change_begin]
  508. X3J13 voted in March 1988 (PUSH-EVALUATION-ORDER)   to clarify order of
  509. evaluation (see section 7.2).
  510. [change_end]
  511.  
  512. [Function]
  513. conjugate number
  514.  
  515. This returns the complex conjugate of number. The conjugate of a non-complex
  516. number is itself. For a complex number z,
  517.  
  518. (conjugate z) == (complex (realpart z) (- (imagpart z)))
  519.  
  520. For example:
  521.  
  522. (conjugate #C(3/5 4/5)) => #C(3/5 -4/5)
  523. (conjugate #C(0.0D0 -1.0D0)) => #C(0.0D0 1.0D0)
  524. (conjugate 3.7) => 3.7
  525.  
  526. [Function]
  527. gcd &rest integers
  528.  
  529. This returns the greatest common divisor of all the arguments, which must be
  530. integers. The result of gcd is always a non-negative integer. If one argument
  531. is given, its absolute value is returned. If no arguments are given, gcd
  532. returns 0, which is an identity for this operation. For three or more
  533. arguments,
  534.  
  535. (gcd a b c ... z) == (gcd (gcd a b) c ... z)
  536.  
  537. Here are some examples of the use of gcd:
  538.  
  539. (gcd 91 -49) => 7
  540. (gcd 63 -42 35) => 7
  541. (gcd 5) => 5
  542. (gcd -4) => 4
  543. (gcd) => 0
  544.  
  545. [Function]
  546. lcm integer &rest more-integers
  547.  
  548. This returns the least common multiple of its arguments, which must be
  549. integers. The result of lcm is always a non-negative integer. For two arguments
  550. that are not both zero,
  551.  
  552. (lcm a b) == (/ (abs (* a b)) (gcd a b))
  553.  
  554. If one or both arguments are zero,
  555.  
  556. (lcm a 0) == (lcm 0 a) == 0
  557.  
  558. For one argument, lcm returns the absolute value of that argument. For three or
  559. more arguments,
  560.  
  561. (lcm a b c ... z) == (lcm (lcm a b) c ... z)
  562.  
  563. Some examples:
  564.  
  565. (lcm 14 35) => 70
  566. (lcm 0 5) => 0
  567. (lcm 1 2 3 4 5 6) => 60
  568.  
  569. Mathematically, (lcm) should return infinity. Because Common Lisp does not have
  570. a representation for infinity, lcm, unlike gcd, always requires at least one
  571. argument.
  572.  
  573. [change_begin]
  574. X3J13 voted in January 1989 (LCM-NO-ARGUMENTS)   to specify that (lcm) => 1.
  575.  
  576. This is one of my biggest boners. The identity for lcm is of course 1, not
  577. infinity, and so (lcm) ought to have been defined to return 1. Sorry about
  578. that, though in point of fact very few users have complained to me that this
  579. mistake in the first edition has cramped their programming style.
  580. [change_end]
  581.  
  582. -------------------------------------------------------------------------------
  583.  
  584. 12.5. Irrational and Transcendental Functions
  585.  
  586. Common Lisp provides no data type that can accurately represent irrational
  587. numerical values. The functions in this section are described as if the results
  588. were mathematically accurate, but actually they all produce floating-point
  589. approximations to the true mathematical result in the general case. In some
  590. places mathematical identities are set forth that are intended to elucidate the
  591. meanings of the functions; however, two mathematically identical expressions
  592. may be computationally different because of errors inherent in the
  593. floating-point approximation process.
  594.  
  595. When the arguments to a function in this section are all rational and the true
  596. mathematical result is also (mathematically) rational, then unless otherwise
  597. noted an implementation is free to return either an accurate result of type
  598. rational or a single-precision floating-point approximation. If the arguments
  599. are all rational but the result cannot be expressed as a rational number, then
  600. a single-precision floating-point approximation is always returned.
  601.  
  602. [change_begin]
  603. X3J13 voted in March 1989 (COMPLEX-RATIONAL-RESULT)   to clarify that the
  604. provisions of the previous paragraph apply to complex numbers. If the arguments
  605. to a function are all of type (or rational (complex rational)) and the true
  606. mathematical result is (mathematically) a complex number with rational real and
  607. imaginary parts, then unless otherwise noted an implementation is free to
  608. return either an accurate result of type (or rational (complex rational)) or a
  609. single-precision floating-point approximation of type single-float (permissible
  610. only if the imaginary part of the true mathematical result is zero) or (complex
  611. single-float). If the arguments are all of type (or rational (complex
  612. rational)) but the result cannot be expressed as a rational or complex rational
  613. number, then the returned value will be of type single-float (permissible only
  614. if the imaginary part of the true mathematical result is zero) or (complex
  615. single-float).
  616. [change_end]
  617.  
  618. The rules of floating-point contagion and complex contagion are effectively
  619. obeyed by all the functions in this section except expt, which treats some
  620. cases of rational exponents specially. When, possibly after contagious
  621. conversion, all of the arguments are of the same floating-point or complex
  622. floating-point type, then the result will be of that same type unless otherwise
  623. noted.
  624.  
  625. -------------------------------------------------------------------------------
  626. Implementation note: There is a ``floating-point cookbook'' by Cody and Waite
  627. [14] that may be a useful aid in implementing the functions defined in this
  628. section.
  629. -------------------------------------------------------------------------------
  630.  
  631. -------------------------------------------------------------------------------
  632.  
  633.    *  Exponential and Logarithmic Functions
  634.    *  Trigonometric and Related Functions
  635.    *  Branch Cuts, Principal Values, and Boundary Conditions in the Complex
  636.      Plane
  637.  
  638. -------------------------------------------------------------------------------
  639.  
  640. 12.5.1. Exponential and Logarithmic Functions
  641.  
  642. Along with the usual one-argument and two-argument exponential and logarithm
  643. functions, sqrt is considered to be an exponential function, because it raises
  644. a number to the power 1/2.
  645.  
  646. [Function]
  647. exp number
  648.  
  649. Returns e raised to the power number, where e is the base of the natural
  650. logarithms.
  651.  
  652. [Function]
  653. expt base-number power-number
  654.  
  655. Returns base-number raised to the power power-number. If the base-number is of
  656. type rational and the power-number is an integer, the calculation will be exact
  657. and the result will be of type rational; otherwise a floating-point
  658. approximation may result.
  659.  
  660. [change_begin]
  661. X3J13 voted in March 1989 (COMPLEX-RATIONAL-RESULT)   to clarify that
  662. provisions similar to those of the previous paragraph apply to complex numbers.
  663. If the base-number is of type (complex rational) and the power-number is an
  664. integer, the calculation will also be exact and the result will be of type (or
  665. rational (complex rational)); otherwise a floating-point or complex
  666. floating-point approximation may result.
  667. [change_end]
  668.  
  669. When power-number is 0 (a zero of type integer), then the result is always the
  670. value 1 in the type of base-number, even if the base-number is zero (of any
  671. type). That is:
  672.  
  673. (expt x 0) == (coerce 1 (type-of x))
  674.  
  675. If the power-number is a zero of any other data type, then the result is also
  676. the value 1, in the type of the arguments after the application of the
  677. contagion rules, with one exception: it is an error if base-number is zero when
  678. the power-number is a zero not of type integer.
  679.  
  680. Implementations of expt are permitted to use different algorithms for the cases
  681. of a rational power-number and a floating-point power-number; the motivation is
  682. that in many cases greater accuracy can be achieved for the case of a rational
  683. power-number. For example, (expt pi 16) and (expt pi 16.0) may yield slightly
  684. different results if the first case is computed by repeated squaring and the
  685. second by the use of logarithms. Similarly, an implementation might choose to
  686. compute (expt x 3/2) as if it had been written (sqrt (expt x 3)), perhaps
  687. producing a more accurate result than would (expt x 1.5). It is left to the
  688. implementor to determine the best strategies.
  689.  
  690. [change_begin]
  691. X3J13 voted in January 1989 (EXPT-RATIO)   to clarify that the preceding remark
  692. is in error, because (sqrt (expt x 3)) does not produce the same value as (expt
  693. x 3/2) in most cases, and to specify that the specification of the principal
  694. value of expt as given in section 12.5.3 should be regarded as definitive.
  695.  
  696. As an example of the difficulty, let   . Then   , but   . Another example is
  697. x=-1; then   , but   .
  698. [change_end]
  699.  
  700. The result of expt can be a complex number, even when neither argument is
  701. complex, if base-number is negative and power-number is not an integer. The
  702. result is always the principal complex value. Note that (expt -8 1/3) is not
  703. permitted to return -2; while -2 is indeed one of the cube roots of -8, it is
  704. not the principal cube root, which is a complex number approximately equal to
  705. #C(1.0 1.73205).
  706.  
  707. [change_begin]
  708. Notice of correction. The first edition gave the incorrect value #C(0.5
  709. 1.73205) for the principal cube root of -8. The correct value is #C(1.0
  710. 1.73205), that is, 1+SQRT(3)i. I simply don't know what I was thinking of!
  711. [change_end]
  712.  
  713. [Function]
  714. log number &optional base
  715.  
  716. Returns the logarithm of number in the base base, which defaults to e, the base
  717. of the natural logarithms. For example:
  718.  
  719. (log 8.0 2) => 3.0
  720. (log 100.0 10) => 2.0
  721.  
  722. The result of (log 8 2) may be either 3 or 3.0, depending on the
  723. implementation.
  724.  
  725. Note that log may return a complex result when given a non-complex argument if
  726. the argument is negative. For example:
  727.  
  728. (log -1.0) == (complex 0.0 (float pi 0.0))
  729.  
  730. [change_begin]
  731. X3J13 voted in January 1989 (IEEE-ATAN-BRANCH-CUT)   to specify certain
  732. floating-point behavior when minus zero is supported. As a part of that vote it
  733. approved a mathematical definition of complex logarithm in terms of real
  734. logarithm, absolute value, arc tangent of two real arguments, and the phase
  735. function as
  736.  
  737.  
  738. Logarithm               log|z| + i phase z
  739.  
  740. This specifies the branch cuts precisely whether minus zero is supported or
  741. not; see phase and atan.
  742. [change_end]
  743.  
  744. [Function]
  745. sqrt number
  746.  
  747. Returns the principal square root of number. If the number is not complex but
  748. is negative, then the result will be a complex number. For example:
  749.  
  750. (sqrt 9.0) => 3.0
  751. (sqrt -9.0) => #c(0.0 3.0)
  752.  
  753. The result of (sqrt 9) may be either 3 or 3.0, depending on the implementation.
  754. The result of (sqrt -9) may be either #c(0 3) or #c(0.0 3.0).
  755.  
  756. [change_begin]
  757. X3J13 voted in January 1989 (IEEE-ATAN-BRANCH-CUT)   to specify certain
  758. floating-point behavior when minus zero is supported. As a part of that vote it
  759. approved a mathematical definition of complex square root in terms of complex
  760. logarithm and exponential functions as
  761.  
  762.  
  763. Square root
  764.  
  765. This specifies the branch cuts precisely whether minus zero is supported or
  766. not; see phase and atan.
  767. [change_end]
  768.  
  769. [Function]
  770. isqrt integer
  771.  
  772. Integer square root: the argument must be a non-negative integer, and the
  773. result is the greatest integer less than or equal to the exact positive square
  774. root of the argument. For example:
  775.  
  776. (isqrt 9) => 3
  777. (isqrt 12) => 3
  778. (isqrt 300) => 17
  779. (isqrt 325) => 18
  780.  
  781. -------------------------------------------------------------------------------
  782.  
  783. 12.5.2. Trigonometric and Related Functions
  784.  
  785. Some of the functions in this section, such as abs and signum, are apparently
  786. unrelated to trigonometric functions when considered as functions of real
  787. numbers only. The way in which they are extended to operate on complex numbers
  788. makes the trigonometric connection clear.
  789.  
  790. [Function]
  791. abs number
  792.  
  793. Returns the absolute value of the argument. For a non-complex number x,
  794.  
  795. (abs x) == (if (minusp x) (- x) x)
  796.  
  797. and the result is always of the same type as the argument.
  798.  
  799. For a complex number z, the absolute value may be computed as
  800.  
  801. (sqrt (+ (expt (realpart z) 2) (expt (imagpart z) 2)))
  802.  
  803. -------------------------------------------------------------------------------
  804. Implementation note: The careful implementor will not use this formula directly
  805. for all complex numbers but will instead handle very large or very small
  806. components specially to avoid intermediate overflow or underflow.
  807. -------------------------------------------------------------------------------
  808.  
  809. For example:
  810.  
  811. (abs #c(3.0 -4.0)) => 5.0
  812.  
  813. The result of (abs #c(3 4)) may be either 5 or 5.0, depending on the
  814. implementation.
  815.  
  816. [Function]
  817. phase number
  818.  
  819. The phase of a number is the angle part of its polar representation as a
  820. complex number. That is,
  821.  
  822. (phase z) == (atan (imagpart z) (realpart z))
  823.  
  824. [old_change_begin]
  825. The result is in radians, in the range -pi (exclusive) to +pi (inclusive). The
  826. phase of a positive non-complex number is zero; that of a negative non-complex
  827. number is pi. The phase of zero is arbitrarily defined to be zero.
  828. [old_change_end]
  829.  
  830. [change_begin]
  831. X3J13 voted in January 1989 (IEEE-ATAN-BRANCH-CUT)   to specify certain
  832. floating-point behavior when minus zero is supported; phase is still defined in
  833. terms of atan as above, but thanks to a change in atan the range of phase
  834. becomes -pi inclusive to +pi inclusive. The value - results from an argument
  835. whose real part is negative and whose imaginary part is minus zero. The phase
  836. function therefore has a branch cut along the negative real axis. The phase of
  837. +0+0i is +0, of +0-0i is -0, of -0+0i is +pi, and of -0-0i is -pi.
  838. [change_end]
  839.  
  840. If the argument is a complex floating-point number, the result is a
  841. floating-point number of the same type as the components of the argument. If
  842. the argument is a floating-point number, the result is a floating-point number
  843. of the same type. If the argument is a rational number or complex rational
  844. number, the result is a single-format floating-point number.
  845.  
  846. [Function]
  847. signum number
  848.  
  849. By definition,
  850.  
  851. (signum x) == (if (zerop x) x (/ x (abs x)))
  852.  
  853. For a rational number, signum will return one of -1, 0, or 1 according to
  854. whether the number is negative, zero, or positive. For a floating-point number,
  855. the result will be a floating-point number of the same format whose value is
  856. -1, 0, or 1. For a complex number z, (signum z) is a complex number of the same
  857. phase but with unit magnitude, unless z is a complex zero, in which case the
  858. result is z. For example:
  859.  
  860. (signum 0) => 0
  861. (signum -3.7L5) => -1.0L0
  862. (signum 4/5) => 1
  863. (signum #C(7.5 10.0)) => #C(0.6 0.8)
  864. (signum #C(0.0 -14.7)) => #C(0.0 -1.0)
  865.  
  866. For non-complex rational numbers, signum is a rational function, but it may be
  867. irrational for complex arguments.
  868.  
  869. [Function]
  870. sin radians
  871. cos radians
  872. tan radians
  873.  
  874. sin returns the sine of the argument, cos the cosine, and tan the tangent. The
  875. argument is in radians. The argument may be complex.
  876.  
  877. [Function]
  878. cis radians
  879.  
  880. This computes   . The name cis means ``cos + i sin,'' because   . The argument
  881. is in radians and may be any non-complex number. The result is a complex number
  882. whose real part is the cosine of the argument and whose imaginary part is the
  883. sine. Put another way, the result is a complex number whose phase is the equal
  884. to the argument (mod 2) and whose magnitude is unity.
  885.  
  886. -------------------------------------------------------------------------------
  887. Implementation note: Often it is cheaper to calculate the sine and cosine of a
  888. single angle together than to perform two disjoint calculations.
  889. -------------------------------------------------------------------------------
  890.  
  891. [Function]
  892. asin number
  893. acos number
  894.  
  895. asin returns the arc sine of the argument, and acos the arc cosine. The result
  896. is in radians. The argument may be complex.
  897.  
  898. The arc sine and arc cosine functions may be defined mathematically for an
  899. argument z as follows:
  900.  
  901.  
  902. Arc sine
  903.  
  904. Arc cosine
  905.  
  906. Note that the result of asin or acos may be complex even if the argument is not
  907. complex; this occurs when the absolute value of the argument is greater than 1.
  908.  
  909. [change_begin]
  910. Kahan [25] suggests for acos the defining formula
  911.  
  912.  
  913. Arc cosine
  914.  
  915. or even the much simpler   . Both equations are mathematically equivalent to
  916. the formula shown above.
  917. [change_end]
  918.  
  919. -------------------------------------------------------------------------------
  920. Implementation note: These formulae are mathematically correct, assuming
  921. completely accurate computation. They may be terrible methods for
  922. floating-point computation. Implementors should consult a good text on
  923. numerical analysis. The formulae given above are not necessarily the simplest
  924. ones for real-valued computations, either; they are chosen to define the branch
  925. cuts in desirable ways for the complex case.
  926. -------------------------------------------------------------------------------
  927.  
  928. [Function]
  929. atan y &optional x
  930.  
  931. An arc tangent is calculated and the result is returned in radians.
  932.  
  933. With two arguments y and x, neither argument may be complex. The result is the
  934. arc tangent of the quantity y/x. The signs of y and x are used to derive
  935. quadrant information; moreover, x may be zero provided y is not zero. The value
  936. of atan is always between -pi (exclusive) and +pi (inclusive). The following
  937. table details various special cases.
  938.  
  939.  
  940.  
  941. [change_begin]
  942. X3J13 voted in January 1989 (IEEE-ATAN-BRANCH-CUT)   to specify certain
  943. floating-point behavior when minus zero is supported. When there is a minus
  944. zero, the preceding table must be modified slightly:
  945.  
  946.  
  947.  
  948. Note that the case y=0,x=0 is an error in the absence of minus zero, but the
  949. four cases y=0,x=0 are defined in the presence of minus zero.
  950. [change_end]
  951.  
  952. [old_change_begin]
  953. With only one argument y, the argument may be complex. The result is the arc
  954. tangent of y, which may be defined by the following formula:
  955.  
  956. Arc tangent
  957.  
  958. [old_change_end]
  959.  
  960. -------------------------------------------------------------------------------
  961. Implementation note: This formula is mathematically correct, assuming
  962. completely accurate computation. It may be a terrible method for floating-point
  963. computation. Implementors should consult a good text on numerical analysis. The
  964. formula given above is not necessarily the simplest one for real-valued
  965. computations, either; it is chosen to define the branch cuts in desirable ways
  966. for the complex case.
  967. -------------------------------------------------------------------------------
  968.  
  969. [change_begin]
  970. X3J13 voted in January 1989 (COMPLEX-ATAN-BRANCH-CUT)   to replace the
  971. preceding formula with the formula
  972.  
  973.                         log(1+iy) - log(1-iy)
  974. Arc tangent             ---------------------
  975.                                   2i
  976.  
  977. This change alters the direction of continuity for the branch cuts, which
  978. alters the result returned by atan only for arguments on the imaginary axis
  979. that are of magnitude greater than 1. See section 12.5.3 for further details.
  980. [change_end]
  981.  
  982. For a non-complex argument y, the result is non-complex and lies between    and
  983.    (both exclusive).
  984.  
  985. -------------------------------------------------------------------------------
  986. Compatibility note: MacLisp has a function called atan whose range is from 0 to
  987. 2. Almost every other programming language (ANSI Fortran, IBM PL/1, Interlisp)
  988. has a two-argument arc tangent function with range -pi to +pi. Lisp Machine
  989. Lisp provides two two-argument arc tangent functions, atan (compatible with
  990. MacLisp) and atan2 (compatible with all others).
  991.  
  992. Common Lisp makes two-argument atan the standard one with range -pi to +pi.
  993. Observe that this makes the one-argument and two-argument versions of atan
  994. compatible in the sense that the branch cuts do not fall in different places.
  995. The Interlisp one-argument function arctan has a range from 0 to pi, while
  996. nearly every other programming language provides the range    to    for
  997. one-argument arc tangent! Nevertheless, since Interlisp uses the standard
  998. two-argument version of arc tangent, its branch cuts are inconsistent anyway.
  999. -------------------------------------------------------------------------------
  1000.  
  1001. [Constant]
  1002. pi
  1003.  
  1004. This global variable has as its value the best possible approximation to pi in
  1005. long floating-point format. For example:
  1006.  
  1007. (defun sind (x)     ;The argument is in degrees
  1008.   (sin (* x (/ (float pi x) 180))))
  1009.  
  1010. An approximation to pi in some other precision can be obtained by writing
  1011. (float pi x), where x is a floating-point number of the desired precision, or
  1012. by writing (coerce pi type), where type is the name of the desired type, such
  1013. as short-float.
  1014.  
  1015. [Function]
  1016. sinh number
  1017. cosh number
  1018. tanh number
  1019. asinh number
  1020. acosh number
  1021. atanh number
  1022.  
  1023. [old_change_begin]
  1024. These functions compute the hyperbolic sine, cosine, tangent, arc sine, arc
  1025. cosine, and arc tangent functions, which are mathematically defined for an
  1026. argument z as follows:
  1027.  
  1028. Hyperbolic sine
  1029. Hyperbolic cosine
  1030. Hyperbolic tangent
  1031. Hyperbolic arc sine
  1032. Hyperbolic arc cosine
  1033. Hyperbolic arc tangent           WRONG!
  1034.  
  1035. [old_change_end]
  1036.  
  1037. [change_begin]
  1038. WARNING! The formula shown above for hyperbolic arc tangent is incorrect. It is
  1039. not a matter of incorrect branch cuts; it simply does not compute anything like
  1040. a hyperbolic arc tangent. This unfortunate error in the first edition was the
  1041. result of mistranscribing a (correct) APL formula from Penfield's paper [36].
  1042. The formula should have been transcribed as
  1043.  
  1044.  
  1045. Hyperbolic arc tangent
  1046.  
  1047. A proposal was submitted to X3J13 in September 1989 to replace the formulae for
  1048. acosh and atanh. See section 12.5.3 for further discussion.
  1049. [change_end]
  1050.  
  1051. Note that the result of acosh may be complex even if the argument is not
  1052. complex; this occurs when the argument is less than 1. Also, the result of
  1053. atanh may be complex even if the argument is not complex; this occurs when the
  1054. absolute value of the argument is greater than 1.
  1055.  
  1056. -------------------------------------------------------------------------------
  1057. Implementation note: These formulae are mathematically correct, assuming
  1058. completely accurate computation. They may be terrible methods for
  1059. floating-point computation. Implementors should consult a good text on
  1060. numerical analysis. The formulae given above are not necessarily the simplest
  1061. ones for real-valued computations, either; they are chosen to define the branch
  1062. cuts in desirable ways for the complex case.
  1063. -------------------------------------------------------------------------------
  1064.  
  1065. -------------------------------------------------------------------------------
  1066.  
  1067. 12.5.3. Branch Cuts, Principal Values, and Boundary Conditions in the Complex
  1068. Plane
  1069.  
  1070. Many of the irrational and transcendental functions are multiply defined in the
  1071. complex domain; for example, there are in general an infinite number of complex
  1072. values for the logarithm function. In each such case, a principal value must be
  1073. chosen for the function to return. In general, such values cannot be chosen so
  1074. as to make the range continuous; lines in the domain called branch cuts must be
  1075. defined, which in turn define the discontinuities in the range.
  1076.  
  1077. Common Lisp defines the branch cuts, principal values, and boundary conditions
  1078. for the complex functions following a proposal for complex functions in APL
  1079. [36]. The contents of this section are borrowed largely from that proposal.
  1080.  
  1081. -------------------------------------------------------------------------------
  1082. Compatibility note: The branch cuts defined here differ in a few very minor
  1083. respects from those advanced by W. Kahan, who considers not only the ``usual''
  1084. definitions but also the special modifications necessary for IEEE proposed
  1085. floating-point arithmetic, which has infinities and minus zero as explicit
  1086. computational objects. For example, he proposes that SQRT(-4+0i)=2i, but
  1087. SQRT(-4-0i)=-2i.
  1088.  
  1089. It may be that the differences between the APL proposal and Kahan's proposal
  1090. will be ironed out. If so, Common Lisp may be changed as necessary to be
  1091. compatible with these other groups. Any changes from the specification below
  1092. are likely to be quite minor, probably concerning primarily questions of which
  1093. side of a branch cut is continuous with the cut itself.
  1094. -------------------------------------------------------------------------------
  1095.  
  1096. [change_begin]
  1097. Indeed, X3J13 voted in January 1989 (COMPLEX-ATAN-BRANCH-CUT)   to alter the
  1098. direction of continuity for the branch cuts of atan, and also
  1099. (IEEE-ATAN-BRANCH-CUT)   to address the treatment of branch cuts in
  1100. implementations that have a distinct floating-point minus zero.
  1101.  
  1102. The treatment of minus zero centers in two-argument atan. If there is no minus
  1103. zero, then the branch cut runs just below the negative real axis as before, and
  1104. the range of two-argument atan is   . If there is a minus zero, however, then
  1105. the branch cut runs precisely on the negative real axis, skittering between
  1106. pairs of numbers of the form -x+/-0i, and the range of two-argument atan is   .
  1107.  
  1108. The treatment of minus zero by all other irrational and transcendental
  1109. functions is then specified by defining those functions in terms of
  1110. two-argument atan. First, phase is defined in terms of two-argument atan, and
  1111. complex abs in terms of real sqrt; then complex log is defined in terms of
  1112. phase, abs, and real log; then complex sqrt in terms of complex log; and
  1113. finally all others are defined in terms of these.
  1114.  
  1115. Kahan [25] treats these matters in some detail and also suggests specific
  1116. algorithms for implementing irrational and transcendental functions in IEEE
  1117. standard floating-point arithmetic [23].
  1118.  
  1119. Remarks in the first edition about the direction of the continuity of branch
  1120. cuts continue to hold in the absence of minus zero and may be ignored if minus
  1121. zero is supported; since all branch cuts happen to run along the principal
  1122. axes, they run between plus zero and minus zero, and so each sort of zero is
  1123. associated with the obvious quadrant.
  1124. [change_end]
  1125.  
  1126. sqrt
  1127.      The branch cut for square root lies along the negative real axis,
  1128.      continuous with quadrant II. The range consists of the right half-plane,
  1129.      including the non-negative imaginary axis and excluding the negative
  1130.      imaginary axis.
  1131.  
  1132. [change_begin]
  1133.  
  1134.      X3J13 voted in January 1989 (IEEE-ATAN-BRANCH-CUT)   to specify certain
  1135.      floating-point behavior when minus zero is supported. As a part of that
  1136.      vote it approved a mathematical definition of complex square root:
  1137.  
  1138.  
  1139.  
  1140.      This defines the branch cuts precisely, whether minus zero is supported or
  1141.      not.
  1142.  
  1143. [change_end]
  1144.  
  1145. phase
  1146.      The branch cut for the phase function lies along the negative real axis,
  1147.      continuous with quadrant II. The range consists of that portion of the
  1148.      real axis between -pi (exclusive) and pi (inclusive).
  1149.  
  1150. [change_begin]
  1151.  
  1152.      X3J13 voted in January 1989 (IEEE-ATAN-BRANCH-CUT)   to specify certain
  1153.      floating-point behavior when minus zero is supported. As a part of that
  1154.      vote it approved a mathematical definition of phase:
  1155.  
  1156.  
  1157.  
  1158.      where Fz is the imaginary part of z and Rz the real part of z. This
  1159.      defines the branch cuts precisely, whether minus zero is supported or not.
  1160.  
  1161. [change_end]
  1162.  
  1163. log  The branch cut for the logarithm function of one argument (natural
  1164.      logarithm) lies along the negative real axis, continuous with quadrant II.
  1165.      The domain excludes the origin. For a complex number z, z is defined to be
  1166.  
  1167.       log z = (log|z|) + i(phase z)
  1168.  
  1169.      Therefore the range of the one-argument logarithm function is that strip
  1170.      of the complex plane containing numbers with imaginary parts between -pi
  1171.      (exclusive) and pi (inclusive).
  1172.  
  1173.      [change_begin]
  1174.      The X3J13 vote on minus zero (IEEE-ATAN-BRANCH-CUT)   would alter that
  1175.      exclusive bound of -pi to be inclusive if minus zero is supported.
  1176.      [change_end]
  1177.  
  1178.      The two-argument logarithm function is defined as   . This defines the
  1179.      principal values precisely. The range of the two-argument logarithm
  1180.      function is the entire complex plane. It is an error if z is zero. If z is
  1181.      non-zero and b is zero, the logarithm is taken to be zero.
  1182.  
  1183. exp  The simple exponential function has no branch cut.
  1184.  
  1185. expt
  1186.      The two-argument exponential function is defined as   . This defines the
  1187.      principal values precisely. The range of the two-argument exponential
  1188.      function is the entire complex plane. Regarded as a function of x, with b
  1189.      fixed, there is no branch cut. Regarded as a function of b, with x fixed,
  1190.      there is in general a branch cut along the negative real axis, continuous
  1191.      with quadrant II. The domain excludes the origin. By definition,   . If
  1192.      b=0 and the real part of x is strictly positive, then   . For all other
  1193.      values of x,    is an error.
  1194.  
  1195. asin
  1196.      The following definition for arc sine determines the range and branch
  1197.      cuts:
  1198.  
  1199.  
  1200.  
  1201.      This is equivalent to the formula
  1202.  
  1203.  
  1204.  
  1205.      recommended by Kahan [25].
  1206.  
  1207.      The branch cut for the arc sine function is in two pieces: one along the
  1208.      negative real axis to the left of -1 (inclusive), continuous with quadrant
  1209.      II, and one along the positive real axis to the right of 1 (inclusive),
  1210.      continuous with quadrant IV. The range is that strip of the complex plane
  1211.      containing numbers whose real part is between    and   . A number with
  1212.      real part equal to    is in the range if and only if its imaginary part is
  1213.      non-negative; a number with real part equal to    is in the range if and
  1214.      only if its imaginary part is non-positive.
  1215.  
  1216. acos
  1217.      The following definition for arc cosine determines the range and branch
  1218.      cuts:
  1219.  
  1220.  
  1221.  
  1222.      or, which is equivalent,
  1223.  
  1224.  
  1225.  
  1226.      The branch cut for the arc cosine function is in two pieces: one along the
  1227.      negative real axis to the left of -1 (inclusive), continuous with quadrant
  1228.      II, and one along the positive real axis to the right of 1 (inclusive),
  1229.      continuous with quadrant IV. This is the same branch cut as for arc sine.
  1230.      The range is that strip of the complex plane containing numbers whose real
  1231.      part is between zero and pi. A number with real part equal to zero is in
  1232.      the range if and only if its imaginary part is non-negative; a number with
  1233.      real part equal to pi is in the range if and only if its imaginary part is
  1234.      non-positive.
  1235.  
  1236. atan
  1237.      The following definition for (one-argument) arc tangent determines the
  1238.      range and branch cuts:
  1239.  
  1240.  
  1241.  
  1242. [old_change_begin]
  1243.  
  1244.      Beware of simplifying this formula; ``obvious'' simplifications are likely
  1245.      to alter the branch cuts or the values on the branch cuts incorrectly.
  1246.  
  1247.      The branch cut for the arc tangent function is in two pieces: one along
  1248.      the positive imaginary axis above i (exclusive), continuous with quadrant
  1249.      II, and one along the negative imaginary axis below -i (exclusive),
  1250.      continuous with quadrant IV. The points i and -i are excluded from the
  1251.      domain. The range is that strip of the complex plane containing numbers
  1252.      whose real part is between    and   . A number with real part equal to
  1253.      is in the range if and only if its imaginary part is strictly positive; a
  1254.      number with real part equal to    is in the range if and only if its
  1255.      imaginary part is strictly negative. Thus the range of the arc tangent
  1256.      function is identical to that of the arc sine function with the points
  1257.      and    excluded.
  1258.  
  1259. [old_change_end]
  1260.  
  1261. [change_begin]
  1262.  
  1263.      X3J13 voted in January 1989 (COMPLEX-ATAN-BRANCH-CUT)   to replace the
  1264.      formula shown above with the formula
  1265.  
  1266.  
  1267.  
  1268.      This is equivalent to the formula
  1269.  
  1270.  
  1271.  
  1272.      recommended by Kahan [25]. It causes the upper branch cut to be continuous
  1273.      with quadrant I rather than quadrant II, and the lower branch cut to be
  1274.      continuous with quadrant III rather than quadrant IV; otherwise it agrees
  1275.      with the formula of the first edition. Therefore this change alters the
  1276.      result returned by atan only for arguments on the positive imaginary axis
  1277.      that are of magnitude greater than 1. The full description for this new
  1278.      formula is as follows.
  1279.  
  1280.      The branch cut for the arc tangent function is in two pieces: one along
  1281.      the positive imaginary axis above i (exclusive), continuous with quadrant
  1282.      I, and one along the negative imaginary axis below -i (exclusive),
  1283.      continuous with quadrant III. The points i and -i are excluded from the
  1284.      domain. The range is that strip of the complex plane containing numbers
  1285.      whose real part is between    and   . A number with real part equal to
  1286.      is in the range if and only if its imaginary part is strictly negative; a
  1287.      number with real part equal to    is in the range if and only if its
  1288.      imaginary part is strictly positive. Thus the range of the arc tangent
  1289.      function is not identical to that of the arc sine function.
  1290.  
  1291. [change_end]
  1292.  
  1293. asinh
  1294.      The following definition for the inverse hyperbolic sine determines the
  1295.      range and branch cuts:
  1296.  
  1297.  
  1298.  
  1299.      The branch cut for the inverse hyperbolic sine function is in two pieces:
  1300.      one along the positive imaginary axis above i (inclusive), continuous with
  1301.      quadrant I, and one along the negative imaginary axis below -i
  1302.      (inclusive), continuous with quadrant III. The range is that strip of the
  1303.      complex plane containing numbers whose imaginary part is between    and
  1304.        . A number with imaginary part equal to    is in the range if and only
  1305.      if its real part is non-positive; a number with imaginary part equal to
  1306.      is in the range if and only if its real part is non-negative.
  1307.  
  1308. acosh
  1309.      The following definition for the inverse hyperbolic cosine determines the
  1310.      range and branch cuts:
  1311.  
  1312.  
  1313.  
  1314. [change_begin]
  1315.  
  1316.      Kahan [25] suggests the formula
  1317.  
  1318.  
  1319.  
  1320.      pointing out that it yields the same principal value but eliminates a
  1321.      gratuitous removable singularity at z=-1. A proposal was submitted to
  1322.      X3J13 in September 1989 to replace the formula acosh with that recommended
  1323.      by Kahan. There is a good possibility that it will be adopted.
  1324.  
  1325. [change_end]
  1326.  
  1327.      The branch cut for the inverse hyperbolic cosine function lies along the
  1328.      real axis to the left of 1 (inclusive), extending indefinitely along the
  1329.      negative real axis, continuous with quadrant II and (between 0 and 1) with
  1330.      quadrant I. The range is that half-strip of the complex plane containing
  1331.      numbers whose real part is non-negative and whose imaginary part is
  1332.      between -pi (exclusive) and pi (inclusive). A number with real part zero
  1333.      is in the range if its imaginary part is between zero (inclusive) and pi
  1334.      (inclusive).
  1335.  
  1336. atanh
  1337.      The following definition for the inverse hyperbolic tangent determines the
  1338.      range and branch cuts:
  1339.  
  1340. [old_change_begin]
  1341.  
  1342.              WRONG!
  1343.  
  1344. [old_change_end]
  1345.  
  1346. [change_begin]
  1347.  
  1348.      WARNING! The formula shown above for hyperbolic arc tangent is incorrect.
  1349.      It is not a matter of incorrect branch cuts; it simply does not compute
  1350.      anything like a hyperbolic arc tangent. This unfortunate error in the
  1351.      first edition was the result of mistranscribing a (correct) APL formula
  1352.      from Penfield's paper [36]. The formula should have been transcribed as
  1353.  
  1354.  
  1355.  
  1356. [change_end]
  1357.  
  1358. [old_change_begin]
  1359.  
  1360.      Beware of simplifying this formula; ``obvious'' simplifications are likely
  1361.      to alter the branch cuts or the values on the branch cuts incorrectly.
  1362.  
  1363.      The branch cut for the inverse hyperbolic tangent function is in two
  1364.      pieces: one along the negative real axis to the left of -1 (inclusive),
  1365.      continuous with quadrant III, and one along the positive real axis to the
  1366.      right of 1 (inclusive), continuous with quadrant I. The points -1 and 1
  1367.      are excluded from the domain. The range is that strip of the complex plane
  1368.      containing numbers whose imaginary part is between    and   . A number
  1369.      with imaginary part equal to    is in the range if and only if its real
  1370.      part is strictly negative; a number with imaginary part equal to    is in
  1371.      the range if and only if its real part is strictly positive. Thus the
  1372.      range of the inverse hyperbolic tangent function is identical to that of
  1373.      the inverse hyperbolic sine function with the points    and    excluded.
  1374.  
  1375. [old_change_end]
  1376.  
  1377. [change_begin]
  1378.  
  1379.      A proposal was submitted to X3J13 in September 1989 to replace the formula
  1380.      atanh with that recommended by Kahan [25]:
  1381.  
  1382.  
  1383.  
  1384.      There is a good possibility that it will be adopted. If it is, the
  1385.      complete description of the branch cuts of atanh will then be as follows.
  1386.  
  1387.      The branch cut for the inverse hyperbolic tangent function is in two
  1388.      pieces: one along the negative real axis to the left of -1 (inclusive),
  1389.      continuous with quadrant II, and one along the positive real axis to the
  1390.      right of 1 (inclusive), continuous with quadrant IV. The points -1 and 1
  1391.      are excluded from the domain. The range is that strip of the complex plane
  1392.      containing numbers whose imaginary part is between    and   . A number
  1393.      with imaginary part equal to    is in the range if and only if its real
  1394.      part is strictly positive; a number with imaginary part equal to    is in
  1395.      the range if and only if its real part is strictly negative. Thus the
  1396.      range of the inverse hyperbolic tangent function is not the same as that
  1397.      of the inverse hyperbolic sine function.
  1398.  
  1399. [change_end]
  1400.  
  1401. With these definitions, the following useful identities are obeyed throughout
  1402. the applicable portion of the complex domain, even on the branch cuts:
  1403.  
  1404. sin iz = i sinh z  sinh iz = i sin z        arctan iz = i arctanh z
  1405. cos iz = cosh z    cosh iz = cos z          arcsinh iz = i arcsin z
  1406. tan iz = i tanh z  arcsin iz = i arcsinh z  arctanh iz = i arctan z
  1407.  
  1408. [change_begin]
  1409. I thought it would be useful to provide some graphs illustrating the behavior
  1410. of the irrational and transcendental functions in the complex plane. It also
  1411. provides an opportunity to show off the Common Lisp code that was used to
  1412. generate them.
  1413.  
  1414. Imagine the complex plane to be decorated as follows. The real and imaginary
  1415. axes are painted with thick lines. Parallels from the axes on both sides at
  1416. distances of 1, 2, and 3 are painted with thin lines; these parallels are
  1417. doubly infinite lines, as are the axes. Four annuli (rings) are painted in
  1418. gradated shades of gray. Ring 1, the inner ring, consists of points whose
  1419. radial distances from the origin lie in the range [1/4, 1/2]; ring 2 is in the
  1420. radial range [3/4, 1]; ring 3, in the range [pi/2, 2]; and ring 4, in the range
  1421. [3, pi]. Ring j is divided into    equal sectors, with each sector painted a
  1422. different shade of gray, darkening as one proceeds counterclockwise from the
  1423. positive real axis.
  1424.  
  1425. We can illustrate the behavior of a numerical function f by considering how it
  1426. maps the complex plane to itself. More specifically, consider each point z of
  1427. the decorated plane. We decorate a new plane by coloring the point f(z) with
  1428. the same color that point z had in the original decorated plane. In other
  1429. words, the newly decorated plane illustrates how the f maps the axes, other
  1430. horizontal and vertical lines, and annuli.
  1431.  
  1432. In each figure we will show only a fragment of the complex plane, with the real
  1433. axis horizontal in the usual manner (-pi to the left, +pi to the right) and the
  1434. imaginary axis vertical (-i below, +pii above). Each fragment shows a region
  1435. containing points whose real and imaginary parts are in the range [-4.1, 4.1].
  1436. The axes of the new plane are shown as very thin lines, with large tick marks
  1437. at integer coordinates and somewhat smaller tick marks at multiples of   .
  1438.  
  1439. Figure 12-1 shows the result of plotting the identity function (quite
  1440. literally); the graph exhibits the decoration of the original plane.
  1441.  
  1442. Figures 12-2 through 12-20 show the graphs for the functions sqrt, exp, log,
  1443. sin, asin, cos, acos, tan, atan, sinh, asinh, cosh, acosh, tanh, and atanh, and
  1444. as a bonus, the graphs for the functions   ,   ,   , and   . All of these are
  1445. related to the trigonometric functions in various ways. For example, if   ,
  1446. then   , and if   , then   . It is instructive to examine the graph for    and
  1447. try to visualize how it transforms the graph for sin into the graph for cos.
  1448.  
  1449. Each figure is accompanied by a commentary on what maps to what and other
  1450. interesting features. None of this material is terribly new; much of it may be
  1451. found in any good textbook on complex analysis. I believe that the particular
  1452. form in which the graphs are presented is novel, as well as the fact that the
  1453. graphs have been generated as PostScript [1] code by Common Lisp code. This
  1454. PostScript code was then fed directly to the typesetting equipment that set the
  1455. pages for this book. Samples of this PostScript code follow the figures
  1456. themselves, after which the code for the entire program is presented.
  1457.  
  1458. In the commentaries that accompany the figures I sometimes speak of mapping the
  1459. points +/- infinity or +/- infinity i. When I say that function f maps
  1460. +infinity to a certain point z, I mean that
  1461.  
  1462.  
  1463.  
  1464. Similarly, when I say that f maps -i to z, I mean that
  1465.  
  1466.  
  1467.  
  1468. In other words, I am considering a limit as one travels out along one of the
  1469. main axes. I also speak in a similar manner of mapping to one of these
  1470. infinities.
  1471.  
  1472. -------------------------------------------------------------------------------
  1473.  
  1474.    *   Identity Plot
  1475.    *   Square Root Function
  1476.    *  Exponential Function
  1477.    *  Natural Logarithm Function
  1478.    *  Function (z-1)/(z+1)
  1479.    *  Function (1+z)/(1-z)
  1480.    *  Sine Function
  1481.    *  Arc Sine Function
  1482.    *  Cosine Function
  1483.    *  Arc Cosine Function
  1484.    *  Tanget Function
  1485.    *  Arc Tangent Function
  1486.    *  Hyperbolic Sine Function
  1487.    *  Hyperbolic Arc Sine Function
  1488.    *  Hyperbolic Cosine Function
  1489.    *  Hyperbolic Arc Cosine Function
  1490.    *  Hyperbolic Tangent Function
  1491.    *   Hyperbolic Arc Tangent Function
  1492.    *  SQRT(1 - SQR(z))
  1493.    *   SQRT(1 + SQR(z))
  1494.  
  1495. Here is a sample of the PostScript code that generated figure 12-1, showing the
  1496. initial scaling, translation, and clipping parameters; the code for one sector
  1497. of the innermost annulus; and the code for the negative imaginary axis. Comment
  1498. lines indicate how path or boundary segments were generated separately and then
  1499. spliced (in order to allow for the places that a singularity might lurk, in
  1500. which case the generating code can ``inch up'' to the problematical argument
  1501. value).
  1502.  
  1503. The size of the entire PostScript file for the identity function was about 68
  1504. kilobytes (2757 lines, including comments). The smallest files were the plots
  1505. for atan and atanh, about 65 kilobytes apiece; the largest were the plots for
  1506. sin, cos, sinh, and cosh, about 138 kilobytes apiece.
  1507.  
  1508. % PostScript file for plot of function IDENTITY
  1509. % Plot is to fit in a region 4.666666666666667 inches square
  1510. %  showing axes extending 4.1 units from the origin.
  1511.  
  1512. 40.97560975609756 40.97560975609756 scale
  1513. 4.1 4.1 translate
  1514. newpath
  1515.   -4.1 -4.1 moveto
  1516.   4.1 -4.1 lineto
  1517.   4.1 4.1 lineto
  1518.   -4.1 4.1 lineto
  1519.   closepath
  1520. clip
  1521. % Moby grid for function IDENTITY
  1522. % Annulus 0.25 0.5 4 0.97 0.45
  1523. % Sector from 4.7124 to 6.2832 (quadrant 3)
  1524. newpath
  1525.   0.0 -0.25 moveto
  1526.   0.0 -0.375 lineto
  1527.   %middle radial
  1528.   0.0 -0.375 lineto
  1529.   0.0 -0.5 lineto
  1530.   %end radial
  1531.   0.0 -0.5 lineto
  1532.   0.092 -0.4915 lineto
  1533.   0.1843 -0.4648 lineto
  1534.   0.273 -0.4189 lineto
  1535.   0.3536 -0.3536 lineto
  1536.   %middle circumferential
  1537.   0.3536 -0.3536 lineto
  1538.   0.413 -0.2818 lineto
  1539.   0.4594 -0.1974 lineto
  1540.   0.4894 -0.1024 lineto
  1541.   0.5 0.0 lineto
  1542.   %end circumferential
  1543.   0.5 0.0 lineto
  1544.   0.375 0.0 lineto
  1545.   %middle radial
  1546.   0.375 0.0 lineto
  1547.   0.25 0.0 lineto
  1548.   %end radial
  1549.   0.25 0.0 lineto
  1550.   0.2297 -0.0987 lineto
  1551.   0.1768 -0.1768 lineto
  1552.   %middle circumferential
  1553.   0.1768 -0.1768 lineto
  1554.   0.0922 -0.2324 lineto
  1555.   0.0 -0.25 lineto
  1556.   %end circumferential
  1557.   closepath
  1558. currentgray   0.45 setgray   fill   setgray
  1559.  
  1560. [2598 lines omitted]
  1561.  
  1562. % Vertical line from (0.0, -0.5) to (0.0, 0.0)
  1563. newpath
  1564.   0.0 -0.5 moveto
  1565.   0.0 0.0 lineto
  1566. 0.05 setlinewidth   1 setlinecap  stroke
  1567. % Vertical line from (0.0, -0.5) to (0.0, -1.0)
  1568. newpath
  1569.   0.0 -0.5 moveto
  1570.   0.0 -1.0 lineto
  1571. 0.05 setlinewidth   1 setlinecap  stroke
  1572. % Vertical line from (0.0, -2.0) to (0.0, -1.0)
  1573. newpath
  1574.   0.0 -2.0 moveto
  1575.   0.0 -1.0 lineto
  1576. 0.05 setlinewidth   1 setlinecap  stroke
  1577. % Vertical line from (0.0, -2.0) to (0.0, -1.1579208923731617E77)
  1578. newpath
  1579.   0.0 -2.0 moveto
  1580.   0.0 -6.3553 lineto
  1581.   0.0 -6.378103166302659 lineto
  1582.   0.0 -6.378103166302659 lineto
  1583.   0.0 -6.378103166302659 lineto
  1584. 0.05 setlinewidth   1 setlinecap  stroke
  1585.  
  1586. [84 lines omitted]
  1587.  
  1588. % End of PostScript file for plot of function IDENTITY
  1589.  
  1590. Here is the program that generated the PostScript code for the graphs shown in
  1591. figures 12-1 through 12-20. It contains a mixture of fairly general mechanisms
  1592. and ad hoc kludges for plotting functions of a single complex argument while
  1593. gracefully handling extremely large and small values, branch cuts,
  1594. singularities, and periodic behavior. The aim was to provide a simple user
  1595. interface that would not require the caller to provide special advice for each
  1596. function to be plotted. The file for figure 12-1, for example, was generated by
  1597. the call (picture 'identity), which resulted in the writing of a file named
  1598. identity-plot.ps.
  1599.  
  1600. The program assumes that any periodic behavior will have a period that is a
  1601. multiple of 2pi; that branch cuts will fall along the real or imaginary axis;
  1602. and that singularities or very large or small values will occur only at the
  1603. origin, at 1 or +/-i, or on the boundaries of the annuli (particularly those
  1604. with radius    or pi). The central function is parametric-path, which accepts
  1605. four arguments: two real numbers that are the endpoints of an interval of real
  1606. numbers, a function that maps this interval into a path in the complex plane,
  1607. and the function to be plotted; the task of parametric-path is to generate
  1608. PostScript code (a series of lineto operations) that will plot an approximation
  1609. to the image of the parametric path as transformed by the function to be
  1610. plotted. Each of the functions hline, vline, -hline, -vline, radial, and
  1611. circumferential takes appropriate parameters and returns a function suitable
  1612. for use as the third argument to parametric-path. There is some code that
  1613. defends against errors (by using ignore-errors) and against certain
  1614. peculiarities of IEEE floating-point arithmetic (the code that checks for
  1615. not-a-number (NaN) results).
  1616.  
  1617. The program is offered here without further comment or apology.
  1618.  
  1619. -------------------------------------------------------------------------------
  1620.  
  1621. (defparameter units-to-show 4.1)
  1622. (defparameter text-width-in-picas 28.0)
  1623. (defparameter device-pixels-per-inch 300)
  1624. (defparameter pixels-per-unit
  1625.   (* (/ (/ text-width-in-picas 6)
  1626.         (* units-to-show 2))
  1627.      device-pixels-per-inch))
  1628.  
  1629. (defparameter big (sqrt (sqrt most-positive-single-float)))
  1630. (defparameter tiny (sqrt (sqrt least-positive-single-float)))
  1631.  
  1632. (defparameter path-really-losing 1000.0)
  1633. (defparameter path-outer-limit (* units-to-show (sqrt 2) 1.1))
  1634. (defparameter path-minimal-delta (/ 10 pixels-per-unit))
  1635. (defparameter path-outer-delta (* path-outer-limit 0.3))
  1636. (defparameter path-relative-closeness 0.00001)
  1637. (defparameter back-off-delta 0.0005)
  1638.  
  1639. (defun comment-line (stream &rest stuff)
  1640.   (format stream "~%% ")
  1641.   (apply #'format stream stuff)
  1642.   (format t "~%% ")
  1643.   (apply #'format t stuff))
  1644.  
  1645. (defun parametric-path (from to paramfn plotfn)
  1646.   (assert (and (plusp from) (plusp to)))
  1647.   (flet ((domainval (x) (funcall paramfn x))
  1648.          (rangeval (x) (funcall plotfn (funcall paramfn x)))
  1649.          (losing (x) (or (null x)
  1650.                          (/= (realpart x) (realpart x))  ;NaN?
  1651.                          (/= (imagpart x) (imagpart x))  ;NaN?
  1652.                          (> (abs (realpart x)) path-really-losing)
  1653.                          (> (abs (imagpart x)) path-really-losing))))
  1654.     (when (> to 1000.0)
  1655.       (let ((f0 (rangeval from))
  1656.             (f1 (rangeval (+ from 1)))
  1657.             (f2 (rangeval (+ from (* 2 pi))))
  1658.             (f3 (rangeval (+ from 1 (* 2 pi))))
  1659.             (f4 (rangeval (+ from (* 4 pi)))))
  1660.         (flet ((close (x y)
  1661.                  (or (< (careful-abs (- x y)) path-minimal-delta)
  1662.                      (< (careful-abs (- x y))
  1663.                         (* (+ (careful-abs x) (careful-abs y))
  1664.                            path-relative-closeness)))))
  1665.           (when (and (close f0 f2)
  1666.                      (close f2 f4)
  1667.                      (close f1 f3)
  1668.                      (or (and (close f0 f1)
  1669.                               (close f2 f3))
  1670.                          (and (not (close f0 f1))
  1671.                               (not (close f2 f3)))))
  1672.             (format t "~&Periodicity detected.")
  1673.             (setq to (+ from (* (signum (- to from)) 2 pi)))))))
  1674.      (let ((fromrange (ignore-errors (rangeval from)))
  1675.           (torange (ignore-errors (rangeval to))))
  1676.       (if (losing fromrange)
  1677.           (if (losing torange)
  1678.               '()
  1679.               (parametric-path (back-off from to) to paramfn plotfn))
  1680.           (if (losing torange)
  1681.               (parametric-path from (back-off to from) paramfn plotfn)
  1682.               (expand-path (refine-path (list from to) #'rangeval)
  1683.                            #'rangeval))))))
  1684.  
  1685. (defun back-off (point other)
  1686.   (if (or (> point 10.0) (< point 0.1))
  1687.       (let ((sp (sqrt point)))
  1688.         (if (or (> point sp other) (< point sp other))
  1689.             sp
  1690.             (* sp (sqrt other))))
  1691.       (+ point (* (signum (- other point)) back-off-delta))))
  1692.  
  1693. (defun careful-abs (z)
  1694.   (cond ((or (> (realpart z) big)
  1695.              (< (realpart z) (- big))
  1696.              (> (imagpart z) big)
  1697.              (< (imagpart z) (- big)))
  1698.          big)
  1699.         ((complexp z) (abs z))
  1700.         ((minusp z) (- z))
  1701.         (t z)))
  1702.  
  1703. (defparameter max-refinements 5000)
  1704.  
  1705. (defun refine-path (original-path rangevalfn)
  1706.   (flet ((rangeval (x) (funcall rangevalfn x)))
  1707.     (let ((path original-path))
  1708.       (do ((j 0 (+ j 1)))
  1709.           ((null (rest path)))
  1710.         (when (zerop (mod (+ j 1) max-refinements))
  1711.               (break "Runaway path"))
  1712.         (let* ((from (first path))
  1713.                (to (second path))
  1714.                (fromrange (rangeval from))
  1715.                (torange (rangeval to))
  1716.                (dist (careful-abs (- torange fromrange)))
  1717.                (mid (* (sqrt from) (sqrt to)))
  1718.                (midrange (rangeval mid)))
  1719.           (cond ((or (and (far-out fromrange) (far-out torange))
  1720.                      (and (< dist path-minimal-delta)
  1721.                           (< (abs (- midrange fromrange))
  1722.                              path-minimal-delta)
  1723.                           ;; Next test is intentionally asymmetric to
  1724.                           ;;  avoid problems with periodic functions.
  1725.                           (< (abs (- (rangeval (/ (+ to (* from 1.5))
  1726.                                                   2.5))
  1727.                                      fromrange))
  1728.                              path-minimal-delta)))
  1729.                  (pop path))
  1730.                 ((= mid from) (pop path))
  1731.                 ((= mid to) (pop path))
  1732.                 (t (setf (rest path) (cons mid (rest path)))))))))
  1733.   original-path)
  1734.  
  1735. (defun expand-path (path rangevalfn)
  1736.   (flet ((rangeval (x) (funcall rangevalfn x)))
  1737.     (let ((final-path (list (rangeval (first path)))))
  1738.       (do ((p (rest path) (cdr p)))
  1739.           ((null p)
  1740.            (unless (rest final-path)
  1741.              (break "Singleton path"))
  1742.            (reverse final-path))
  1743.         (let ((v (rangeval (car p))))
  1744.           (cond ((and (rest final-path)
  1745.                       (not (far-out v))
  1746.                       (not (far-out (first final-path)))
  1747.                       (between v (first final-path)
  1748.                                  (second final-path)))
  1749.                  (setf (first final-path) v))
  1750.                 ((null (rest p))   ;Mustn't omit last point
  1751.                  (push v final-path))
  1752.                 ((< (abs (- v (first final-path))) path-minimal-delta))
  1753.                 ((far-out v)
  1754.                  (unless (and (far-out (first final-path))
  1755.                               (< (abs (- v (first final-path)))
  1756.                                  path-outer-delta))
  1757.                    (push (* 1.01 path-outer-limit (signum v))
  1758.                          final-path)))
  1759.                 (t (push v final-path))))))))
  1760.  
  1761. (defun far-out (x)
  1762.   (> (careful-abs x) path-outer-limit))
  1763.  
  1764. (defparameter between-tolerance 0.000001)
  1765.  
  1766. (defun between (p q r)
  1767.   (let ((px (realpart p)) (py (imagpart p))
  1768.         (qx (realpart q)) (qy (imagpart q))
  1769.         (rx (realpart r)) (ry (imagpart r)))
  1770.     (and (or (<= px qx rx) (>= px qx rx))
  1771.          (or (<= py qy ry) (>= py qy ry))
  1772.          (< (abs (- (* (- qx px) (- ry qy))
  1773.                     (* (- rx qx) (- qy py))))
  1774.             between-tolerance))))
  1775.  
  1776. (defun circle (radius)
  1777.   #'(lambda (angle) (* radius (cis angle))))
  1778.  
  1779. (defun hline (imag)
  1780.   #'(lambda (real) (complex real imag)))
  1781.  
  1782. (defun vline (real)
  1783.   #'(lambda (imag) (complex real imag)))
  1784.  
  1785. (defun -hline (imag)
  1786.   #'(lambda (real) (complex (- real) imag)))
  1787.  
  1788. (defun -vline (real)
  1789.   #'(lambda (imag) (complex real (- imag))))
  1790.  
  1791. (defun radial (phi quadrant)
  1792.   #'(lambda (rho) (repair-quadrant (* rho (cis phi)) quadrant)))
  1793.  
  1794. (defun circumferential (rho quadrant)
  1795.   #'(lambda (phi) (repair-quadrant (* rho (cis phi)) quadrant)))
  1796.  
  1797. ;;; Quadrant is 0, 1, 2, or 3, meaning I, II, III, or IV.
  1798.  
  1799. (defun repair-quadrant (z quadrant)
  1800.   (complex (* (+ (abs (realpart z)) tiny)
  1801.               (case quadrant (0 1.0) (1 -1.0) (2 -1.0) (3 1.0)))
  1802.            (* (+ (abs (imagpart z)) tiny)
  1803.               (case quadrant (0 1.0) (1 1.0) (2 -1.0) (3 -1.0)))))
  1804.  
  1805. (defun clamp-real (x)
  1806.   (if (far-out x)
  1807.       (* (signum x) path-outer-limit)
  1808.       (round-real x)))
  1809.  
  1810. (defun round-real (x)
  1811.   (/ (round (* x 10000.0)) 10000.0))
  1812.  
  1813. (defun round-point (z)
  1814.   (complex (round-real (realpart z)) (round-real (imagpart z))))
  1815.  
  1816. (defparameter hiringshade 0.97)
  1817. (defparameter loringshade 0.45)
  1818.  
  1819. (defparameter ticklength 0.12)
  1820. (defparameter smallticklength 0.09)
  1821.  
  1822. ;;; This determines the pattern of lines and annuli to be drawn.
  1823. (defun moby-grid (&optional (fn 'sqrt) (stream t))
  1824.   (comment-line stream "Moby grid for function ~S" fn)
  1825.   (shaded-annulus 0.25 0.5 4 hiringshade loringshade fn stream)
  1826.   (shaded-annulus 0.75 1.0 8 hiringshade loringshade fn stream)
  1827.   (shaded-annulus (/ pi 2) 2.0 16 hiringshade loringshade fn stream)
  1828.   (shaded-annulus 3 pi 32 hiringshade loringshade fn stream)
  1829.   (moby-lines :horizontal 1.0 fn stream)
  1830.   (moby-lines :horizontal -1.0 fn stream)
  1831.   (moby-lines :vertical 1.0 fn stream)
  1832.   (moby-lines :vertical -1.0 fn stream)
  1833.   (let ((tickline 0.015)
  1834.         (axisline 0.008))
  1835.     (flet ((tick (n) (straight-line (complex n ticklength)
  1836.                                     (complex n (- ticklength))
  1837.                                     tickline
  1838.                                     stream))
  1839.            (smalltick (n) (straight-line (complex n smallticklength)
  1840.                                          (complex n (- smallticklength))
  1841.                                          tickline
  1842.                                          stream)))
  1843.       (comment-line stream "Real axis")
  1844.       (straight-line #c(-5 0) #c(5 0) axisline stream)
  1845.       (dotimes (j (floor units-to-show))
  1846.         (let ((q (+ j 1))) (tick q) (tick (- q))))
  1847.       (dotimes (j (floor units-to-show (/ pi 2)))
  1848.         (let ((q (* (/ pi 2) (+ j 1))))
  1849.           (smalltick q)
  1850.           (smalltick (- q)))))
  1851.     (flet ((tick (n) (straight-line (complex ticklength n)
  1852.                                     (complex (- ticklength) n)
  1853.                                     tickline
  1854.                                     stream))
  1855.            (smalltick (n) (straight-line (complex smallticklength n)
  1856.                                          (complex (- smallticklength) n)
  1857.                                          tickline
  1858.                                          stream)))
  1859.       (comment-line stream "Imaginary axis")
  1860.       (straight-line #c(0 -5) #c(0 5) axisline stream)
  1861.       (dotimes (j (floor units-to-show))
  1862.         (let ((q (+ j 1))) (tick q) (tick (- q))))
  1863.       (dotimes (j (floor units-to-show (/ pi 2)))
  1864.         (let ((q (* (/ pi 2) (+ j 1))))
  1865.           (smalltick q)
  1866.           (smalltick (- q)))))))
  1867.  
  1868. (defun straight-line (from to wid stream)
  1869.   (format stream
  1870.           "~%newpath  ~S ~S moveto  ~S ~S lineto  ~S ~
  1871.            setlinewidth  1  setlinecap  stroke"
  1872.           (realpart from)
  1873.           (imagpart from)
  1874.           (realpart to)
  1875.           (imagpart to)
  1876.  
  1877. ;;; This function draws the lines for the pattern.
  1878. (defun moby-lines (orientation signum plotfn stream)
  1879.   (let ((paramfn (ecase orientation
  1880.                    (:horizontal (if (< signum 0) #'-hline #'hline))
  1881.                    (:vertical (if (< signum 0) #'-vline #'vline)))))
  1882.     (flet ((foo (from to other wid)
  1883.              (ecase orientation
  1884.                (:horizontal
  1885.                 (comment-line stream
  1886.                               "Horizontal line from (~S, ~S) to (~S, ~S)"
  1887.                               (round-real (* signum from))
  1888.                               (round-real other)
  1889.                               (round-real (* signum to))
  1890.                               (round-real other)))
  1891.                (:vertical
  1892.                 (comment-line stream
  1893.                               "Vertical line from (~S, ~S) to (~S, ~S)"
  1894.                               (round-real other)
  1895.                               (round-real (* signum from))
  1896.                               (round-real other)
  1897.                               (round-real (* signum to)))))
  1898.              (postscript-path
  1899.                stream
  1900.                (parametric-path from
  1901.                                 to
  1902.                                 (funcall paramfn other)
  1903.                                 plotfn))
  1904.              (postscript-penstroke stream wid)))
  1905.       (let* ((thick 0.05)
  1906.              (thin 0.02))
  1907.         ;; Main axis
  1908.         (foo 0.5 tiny 0.0 thick)
  1909.         (foo 0.5 1.0 0.0 thick)
  1910.         (foo 2.0 1.0 0.0 thick)
  1911.         (foo 2.0 big 0.0 thick)
  1912.         ;; Parallels at 1 and -1
  1913.         (foo 2.0 tiny 1.0 thin)
  1914.         (foo 2.0 big 1.0 thin)
  1915.         (foo 2.0 tiny -1.0 thin)
  1916.         (foo 2.0 big -1.0 thin)
  1917.         ;; Parallels at 2, 3, -2, -3
  1918.         (foo tiny big 2.0 thin)
  1919.         (foo tiny big -2.0 thin)
  1920.         (foo tiny big 3.0 thin)
  1921.         (foo tiny big -3.0 thin)))))
  1922.  
  1923. (defun splice (p q)
  1924.   (let ((v (car (last p)))
  1925.         (w (first q)))
  1926.     (and (far-out v)
  1927.          (far-out w)
  1928.          (>= (abs (- v w)) path-outer-delta)
  1929.          ;; Two far-apart far-out points.  Try to walk around
  1930.          ;;  outside the perimeter, in the shorter direction.
  1931.          (let* ((pdiff (phase (/ v w)))
  1932.                 (npoints (floor (abs pdiff) (asin .2)))
  1933.                 (delta (/ pdiff (+ npoints 1)))
  1934.                 (incr (cis delta)))
  1935.            (do ((j 0 (+ j 1))
  1936.                 (p (list w "end splice") (cons (* (car p) incr) p)))
  1937.  
  1938. ;;; This function draws the annuli for the pattern.
  1939. (defun shaded-annulus (inner outer sectors firstshade lastshade fn stream)
  1940.   (assert (zerop (mod sectors 4)))
  1941.   (comment-line stream "Annulus ~S ~S ~S ~S ~S"
  1942.                 (round-real inner) (round-real outer)
  1943.                 sectors firstshade lastshade)
  1944.   (dotimes (jj sectors)
  1945.     (let ((j (- sectors jj 1)))
  1946.       (let* ((lophase (+ tiny (* 2 pi (/ j sectors))))
  1947.              (hiphase (* 2 pi (/ (+ j 1) sectors)))
  1948.              (midphase (/ (+ lophase hiphase) 2.0))
  1949.              (midradius (/ (+ inner outer) 2.0))
  1950.              (quadrant (floor (* j 4) sectors)))
  1951.         (comment-line stream "Sector from ~S to ~S (quadrant ~S)"
  1952.                       (round-real lophase)
  1953.                       (round-real hiphase)
  1954.                       quadrant)
  1955.         (let ((p0 (reverse (parametric-path midradius
  1956.                                             inner
  1957.                                             (radial lophase quadrant)
  1958.                                             fn)))
  1959.               (p1 (parametric-path midradius
  1960.                                    outer
  1961.                                    (radial lophase quadrant)
  1962.                                    fn))
  1963.               (p2 (reverse (parametric-path midphase
  1964.                                             lophase
  1965.                                             (circumferential outer
  1966.                                                              quadrant)
  1967.                                             fn)))
  1968.               (p3 (parametric-path midphase
  1969.                                    hiphase
  1970.                                    (circumferential outer quadrant)
  1971.                                    fn))
  1972.               (p4 (reverse (parametric-path midradius
  1973.                                             outer
  1974.                                             (radial hiphase quadrant)
  1975.                                             fn)))
  1976.               (p5 (parametric-path midradius
  1977.                                    inner
  1978.                                    (radial hiphase quadrant)
  1979.                                    fn))
  1980.               (p6 (reverse (parametric-path midphase
  1981.                                             hiphase
  1982.                                             (circumferential inner
  1983.                                                              quadrant)
  1984.                                             fn)))
  1985.               (p7 (parametric-path midphase
  1986.                                    lophase
  1987.                                    (circumferential inner quadrant)
  1988.                                    fn)))
  1989.           (postscript-closed-path stream
  1990.             (append
  1991.               p0 (splice p0 p1) '("middle radial")
  1992.               p1 (splice p1 p2) '("end radial")
  1993.               p2 (splice p2 p3) '("middle circumferential")
  1994.               p3 (splice p3 p4) '("end circumferential")
  1995.               p4 (splice p4 p5) '("middle radial")
  1996.               p5 (splice p5 p6) '("end radial")
  1997.               p6 (splice p6 p7) '("middle circumferential")
  1998.               p7 (splice p7 p0) '("end circumferential")
  1999.               )))
  2000.         (postscript-shade stream
  2001.                           (/ (+ (* firstshade (- (- sectors 1) j))
  2002.                                 (* lastshade j))
  2003.                              (- sectors 1)))))))
  2004.  
  2005. (defun postscript-penstroke (stream wid)
  2006.   (format stream "~%~S setlinewidth   1 setlinecap  stroke"
  2007.           wid))
  2008.  
  2009. (defun postscript-shade (stream shade)
  2010.   (format stream "~%currentgray   ~S setgray   fill   setgray"
  2011.           shade))
  2012.  
  2013. (defun postscript-closed-path (stream path)
  2014.   (unless (every #'far-out (remove-if-not #'numberp path))
  2015.     (postscript-raw-path stream path)
  2016.     (format stream "~%  closepath")))
  2017.  
  2018. (defun postscript-path (stream path)
  2019.   (unless (every #'far-out (remove-if-not #'numberp path))
  2020.     (postscript-raw-path stream path)))
  2021.  
  2022. ;;; Print a path as a series of PostScript "lineto" commands.
  2023. (defun postscript-raw-path (stream path)
  2024.   (format stream "~%newpath")
  2025.   (let ((fmt "~%  ~S ~S moveto"))
  2026.     (dolist (pt path)
  2027.       (cond ((stringp pt)
  2028.              (format stream "~%  %~A" pt))
  2029.             (t (format stream
  2030.                        fmt
  2031.                        (clamp-real (realpart pt))
  2032.                        (clamp-real (imagpart pt)))
  2033.                (setq fmt "~%  ~S ~S lineto"))))))
  2034.  
  2035. ;;; Definitions of functions to be plotted that are not
  2036. ;;; standard Common Lisp functions.
  2037.  
  2038. (defun one-plus-over-one-minus (x) (/ (+ 1 x) (- 1 x)))
  2039.  
  2040. (defun one-minus-over-one-plus (x) (/ (- 1 x) (+ 1 x)))
  2041.  
  2042. (defun sqrt-square-minus-one (x) (sqrt (- 1 (* x x))))
  2043.  
  2044. (defun sqrt-one-plus-square (x) (sqrt (+ 1 (* x x))))
  2045.  
  2046. ;;; Because X3J13 voted for a new definition of the atan function,
  2047. ;;; the following definition was used in place of the atan function
  2048. ;;; provided by the Common Lisp implementation I was using.
  2049.  
  2050. (defun good-atan (x)
  2051.   (/ (- (log (+ 1 (* x #c(0 1))))
  2052.         (log (- 1 (* x #c(0 1)))))
  2053.      #c(0 2)))
  2054.  
  2055. ;;; Because the first edition had an erroneous definition of atanh,
  2056. ;;; the following definition was used in place of the atanh function
  2057. ;;; provided by the Common Lisp implementation I was using.
  2058.  
  2059. (defun really-good-atanh (x)
  2060.   (/ (- (log (+ 1 x))
  2061.         (log (- 1 x)))
  2062.  
  2063. ;;; This is the main procedure that is intended to be called by a user.
  2064. (defun picture (&optional (fn #'sqrt))
  2065.   (with-open-file (stream (concatenate 'string
  2066.                                        (string-downcase (string fn))
  2067.                                        "-plot.ps")
  2068.                           :direction :output)
  2069.     (format stream "% PostScript file for plot of function ~S~%" fn)
  2070.     (format stream "% Plot is to fit in a region ~S inches square~%"
  2071.             (/ text-width-in-picas 6.0))
  2072.     (format stream
  2073.             "%  showing axes extending ~S units from the origin.~%"
  2074.             units-to-show)
  2075.     (let ((scaling (/ (* text-width-in-picas 12) (* units-to-show 2))))
  2076.       (format stream "~%~S ~:*~S scale" scaling))
  2077.     (format stream "~%~S ~:*~S translate" units-to-show)
  2078.     (format stream "~%newpath")
  2079.     (format stream "~%  ~S ~S moveto" (- units-to-show) (- units-to-show))
  2080.     (format stream "~%  ~S ~S lineto" units-to-show (- units-to-show))
  2081.     (format stream "~%  ~S ~S lineto" units-to-show units-to-show)
  2082.     (format stream "~%  ~S ~S lineto" (- units-to-show) units-to-show)
  2083.     (format stream "~%  closepath")
  2084.     (format stream "~%clip")
  2085.     (moby-grid fn stream)
  2086.     (format stream
  2087.             "~%% End of PostScript file for plot of function ~S"
  2088.             fn)
  2089.     (terpri stream)))
  2090.  
  2091. -------------------------------------------------------------------------------
  2092.  
  2093. 12.6. Type Conversions and Component Extractions on Numbers
  2094.  
  2095. While most arithmetic functions will operate on any kind of number, coercing
  2096. types if necessary, the following functions are provided to allow specific
  2097. conversions of data types to be forced when desired.
  2098.  
  2099. [Function]
  2100. float number &optional other
  2101.  
  2102. This converts any non-complex number to a floating-point number. With no second
  2103. argument, if number is already a floating-point number, then number is
  2104. returned; otherwise a single-float is produced. If the argument other is
  2105. provided, then it must be a floating-point number, and number is converted to
  2106. the same format as other. See also coerce.
  2107.  
  2108. [Function]
  2109. rational number
  2110. rationalize number
  2111.  
  2112. Each of these functions converts any non-complex number to a rational number.
  2113. If the argument is already rational, it is returned. The two functions differ
  2114. in their treatment of floating-point numbers.
  2115.  
  2116. rational assumes that the floating-point number is completely accurate and
  2117. returns a rational number mathematically equal to the precise value of the
  2118. floating-point number.
  2119.  
  2120. rationalize assumes that the floating-point number is accurate only to the
  2121. precision of the floating-point representation and may return any rational
  2122. number for which the floating-point number is the best available approximation
  2123. of its format; in doing this it attempts to keep both numerator and denominator
  2124. small.
  2125.  
  2126. It is always the case that
  2127.  
  2128. (float (rational x) x) == x
  2129.  
  2130. and
  2131.  
  2132. (float (rationalize x) x) == x
  2133.  
  2134. That is, rationalizing a floating-point number by either method and then
  2135. converting it back to a floating-point number of the same format produces the
  2136. original number. What distinguishes the two functions is that rational
  2137. typically has a simple, inexpensive implementation, whereas rationalize goes to
  2138. more trouble to produce a result that is more pleasant to view and simpler to
  2139. compute with for some purposes.
  2140.  
  2141. [Function]
  2142. numerator rational
  2143. denominator rational
  2144.  
  2145. These functions take a rational number (an integer or ratio) and return as an
  2146. integer the numerator or denominator of the canonical reduced form of the
  2147. rational. The numerator of an integer is that integer; the denominator of an
  2148. integer is 1. Note that
  2149.  
  2150. (gcd (numerator x) (denominator x)) => 1
  2151.  
  2152. The denominator will always be a strictly positive integer; the numerator may
  2153. be any integer. For example:
  2154.  
  2155. (numerator (/ 8 -6)) => -4
  2156. (denominator (/ 8 -6)) => 3
  2157.  
  2158. There is no fix function in Common Lisp because there are several interesting
  2159. ways to convert non-integral values to integers. These are provided by the
  2160. functions below, which perform not only type conversion but also some
  2161. non-trivial calculations as well.
  2162.  
  2163. [Function]
  2164. floor number &optional divisor
  2165. ceiling number &optional divisor
  2166. truncate number &optional divisor
  2167. round number &optional divisor
  2168.  
  2169. In the simple one-argument case, each of these functions converts its argument
  2170. number (which must not be complex) to an integer. If the argument is already an
  2171. integer, it is returned directly. If the argument is a ratio or floating-point
  2172. number, the functions use different algorithms for the conversion.
  2173.  
  2174. floor converts its argument by truncating toward negative infinity; that is,
  2175. the result is the largest integer that is not larger than the argument.
  2176.  
  2177. ceiling converts its argument by truncating toward positive infinity; that is,
  2178. the result is the smallest integer that is not smaller than the argument.
  2179.  
  2180. truncate converts its argument by truncating toward zero; that is, the result
  2181. is the integer of the same sign as the argument and which has the greatest
  2182. integral magnitude not greater than that of the argument.
  2183.  
  2184. round converts its argument by rounding to the nearest integer; if number is
  2185. exactly halfway between two integers (that is, of the form integer+0.5), then
  2186. it is rounded to the one that is even (divisible by 2).
  2187.  
  2188. The following table shows what the four functions produce when given various
  2189. arguments.
  2190.  
  2191. Argument    floor       ceiling     truncate    round
  2192. ----------------------------------------------------------
  2193.  2.6          2           3           2           3
  2194.  2.5          2           3           2           2
  2195.  2.4          2           3           2           2
  2196.  0.7          0           1           0           1
  2197.  0.3          0           1           0           0
  2198. -0.3         -1           0           0           0
  2199. -0.7         -1           0           0          -1
  2200. -2.4         -3          -2          -2          -2
  2201. -2.5         -3          -2          -2          -2
  2202. -2.6         -3          -2          -2          -3
  2203. ----------------------------------------------------------
  2204.  
  2205. If a second argument divisor is supplied, then the result is the appropriate
  2206. type of rounding or truncation applied to the result of dividing the number by
  2207. the divisor. For example, (floor 5 2) == (floor (/ 5 2)) but is potentially
  2208. more efficient.
  2209.  
  2210. [change_begin]
  2211. This statement is not entirely accurate; one should instead say that (values
  2212. (floor 5 2)) == (values (floor (/ 5 2))), because there is a second value to
  2213. consider, as discussed below. In other words, the first values returned by the
  2214. two forms will be the same, but in general the second values will differ.
  2215. Indeed, we have
  2216.  
  2217. (floor 5 2) => 2 and 1
  2218. (floor (/ 5 2)) => 2 and 1/2
  2219.  
  2220. for this example.
  2221. [change_end]
  2222.  
  2223. The divisor may be any non-complex number.
  2224.  
  2225. [change_begin]
  2226. It is generally accepted that it is an error for the divisor to be zero.
  2227. [change_end]
  2228.  
  2229. The one-argument case is exactly like the two-argument case where the second
  2230. argument is 1.
  2231.  
  2232. [change_begin]
  2233. In other words, the one-argument case returns an integer and fractional part
  2234. for the number: (truncate 5.3) => 5.0 and 0.3, for example.
  2235. [change_end]
  2236.  
  2237. Each of the functions actually returns two values, whether given one or two
  2238. arguments. The second result is the remainder and may be obtained using
  2239. multiple-value-bind and related constructs. If any of these functions is given
  2240. two arguments x and y and produces results q and r, then q y+r=x. The first
  2241. result q is always an integer. The remainder r is an integer if both arguments
  2242. are integers, is rational if both arguments are rational, and is floating-point
  2243. if either argument is floating-point. One consequence is that in the
  2244. one-argument case the remainder is always a number of the same type as the
  2245. argument.
  2246.  
  2247. When only one argument is given, the two results are exact; the mathematical
  2248. sum of the two results is always equal to the mathematical value of the
  2249. argument.
  2250.  
  2251. -------------------------------------------------------------------------------
  2252. Compatibility note: The names of the functions floor, ceiling, truncate, and
  2253. round are more accurate than names like fix that have heretofore been used in
  2254. various Lisp systems. The names used here are compatible with standard
  2255. mathematical terminology (and with PL/1, as it happens). In Fortran ifix means
  2256. truncate. Algol 68 provides round and uses entier to mean floor. In MacLisp,
  2257. fix and ifix both mean floor (one is generic, the other flonum-in/fixnum-out).
  2258. In Interlisp, fix means truncate. In Lisp Machine Lisp, fix means floor and
  2259. fixr means round. Standard Lisp provides a fix function but does not specify
  2260. precisely what it does. The existing usage of the name fix is so confused that
  2261. it seemed best to avoid it altogether.
  2262.  
  2263. The names and definitions given here have recently been adopted by Lisp Machine
  2264. Lisp, and MacLisp and NIL (New Implementation of Lisp) seem likely to follow
  2265. suit.
  2266. -------------------------------------------------------------------------------
  2267.  
  2268. [Function]
  2269. mod number divisor
  2270. rem number divisor
  2271.  
  2272. mod performs the operation floor on its two arguments and returns the second
  2273. result of floor as its only result. Similarly, rem performs the operation
  2274. truncate on its arguments and returns the second result of truncate as its only
  2275. result.
  2276.  
  2277. mod and rem are therefore the usual modulus and remainder functions when
  2278. applied to two integer arguments. In general, however, the arguments may be
  2279. integers or floating-point numbers.
  2280.  
  2281. (mod 13 4) => 1                 (rem 13 4) => 1
  2282. (mod -13 4) => 3                (rem -13 4) => -1
  2283. (mod 13 -4) => -3               (rem 13 -4) => 1
  2284. (mod -13 -4) => -1              (rem -13 -4) => -1
  2285. (mod 13.4 1) => 0.4             (rem 13.4 1) => 0.4
  2286. (mod -13.4 1) => 0.6            (rem -13.4 1) => -0.4
  2287.  
  2288. -------------------------------------------------------------------------------
  2289. Compatibility note: The Interlisp function remainder is essentially equivalent
  2290. to the Common Lisp function rem. The MacLisp function remainder is like rem but
  2291. accepts only integer arguments.
  2292. -------------------------------------------------------------------------------
  2293.  
  2294. [Function]
  2295. ffloor number &optional divisor
  2296. fceiling number &optional divisor
  2297. ftruncate number &optional divisor
  2298. fround number &optional divisor
  2299.  
  2300. These functions are just like floor, ceiling, truncate, and round, except that
  2301. the result (the first result of two) is always a floating-point number rather
  2302. than an integer. It is roughly as if ffloor gave its arguments to floor, and
  2303. then applied float to the first result before passing them both back. In
  2304. practice, however, ffloor may be implemented much more efficiently. Similar
  2305. remarks apply to the other three functions. If the first argument is a
  2306. floating-point number, and the second argument is not a floating-point number
  2307. of longer format, then the first result will be a floating-point number of the
  2308. same type as the first argument. For example:
  2309.  
  2310. (ffloor -4.7) => -5.0 and 0.3
  2311. (ffloor 3.5d0) => 3.0d0 and 0.5d0
  2312.  
  2313. [Function]
  2314. decode-float float
  2315. scale-float float integer
  2316. float-radix float
  2317. float-sign float1 &optional float2
  2318. float-digits float
  2319. float-precision float
  2320. integer-decode-float float
  2321.  
  2322. The function decode-float takes a floating-point number and returns three
  2323. values.
  2324.  
  2325. The first value is a new floating-point number of the same format representing
  2326. the significand; the second value is an integer representing the exponent; and
  2327. the third value is a floating-point number of the same format indicating the
  2328. sign (-1.0 or 1.0). Let b be the radix for the floating-point representation;
  2329. then decode-float divides the argument by an integral power of b so as to bring
  2330. its value between 1/b (inclusive) and 1 (exclusive) and returns the quotient as
  2331. the first value. If the argument is zero, however, the result is equal to the
  2332. absolute value of the argument (that is, if there is a negative zero, its
  2333. significand is considered to be a positive zero).
  2334.  
  2335. The second value of decode-float is the integer exponent e to which b must be
  2336. raised to produce the appropriate power for the division. If the argument is
  2337. zero, any integer value may be returned, provided that the identity shown below
  2338. for scale-float holds.
  2339.  
  2340. The third value of decode-float is a floating-point number, of the same format
  2341. as the argument, whose absolute value is 1 and whose sign matches that of the
  2342. argument.
  2343.  
  2344. The function scale-float takes a floating-point number f (not necessarily
  2345. between 1/b and 1) and an integer k, and returns (* f (expt (float b f) k)).
  2346. (The use of scale-float may be much more efficient than using exponentiation
  2347. and multiplication and avoids intermediate overflow and underflow if the final
  2348. result is representable.)
  2349.  
  2350. Note that
  2351.  
  2352. (multiple-value-bind (signif expon sign)
  2353.                      (decode-float f)
  2354.   (scale-float signif expon))
  2355. == (abs f)
  2356.  
  2357. and
  2358.  
  2359. (multiple-value-bind (signif expon sign)
  2360.                      (decode-float f)
  2361.   (* (scale-float signif expon) sign))
  2362. == f
  2363.  
  2364. The function float-radix returns (as an integer) the radix b of the
  2365. floating-point argument.
  2366.  
  2367. The function float-sign returns a floating-point number z such that z and
  2368. float1 have the same sign and also such that z and float2 have the same
  2369. absolute value. The argument float2 defaults to the value of (float 1 float1);
  2370. (float-sign x) therefore always produces a 1.0 or -1.0 of appropriate format
  2371. according to the sign of x. (Note that if an implementation has distinct
  2372. representations for negative zero and positive zero, then (float-sign -0.0) =>
  2373. -1.0.)
  2374.  
  2375. The function float-digits returns, as a non-negative integer, the number of
  2376. radix-b digits used in the representation of its argument (including any
  2377. implicit digits, such as a ``hidden bit''). The function float-precision
  2378. returns, as a non-negative integer, the number of significant radix-b digits
  2379. present in the argument; if the argument is (a floating-point) zero, then the
  2380. result is (an integer) zero. For normalized floating-point numbers, the results
  2381. of float-digits and float-precision will be the same, but the precision will be
  2382. less than the number of representation digits for a denormalized or zero
  2383. number.
  2384.  
  2385. The function integer-decode-float is similar to decode-float but for its first
  2386. value returns, as an integer, the significand scaled so as to be an integer.
  2387. For an argument f, this integer will be strictly less than
  2388.  
  2389. (expt b (float-precision f))
  2390.  
  2391. but no less than
  2392.  
  2393. (expt b (- (float-precision f) 1))
  2394.  
  2395. except that if f is zero, then the integer value will be zero.
  2396.  
  2397. The second value bears the same relationship to the first value as for
  2398. decode-float:
  2399.  
  2400. (multiple-value-bind (signif expon sign)
  2401.                      (integer-decode-float f)
  2402.   (scale-float (float signif f) expon))
  2403. == (abs f)
  2404.  
  2405. The third value of integer-decode-float will be 1 or -1.
  2406.  
  2407. -------------------------------------------------------------------------------
  2408. Rationale: These functions allow the writing of machine-independent, or at
  2409. least machine-parameterized, floating-point software of reasonable efficiency.
  2410. -------------------------------------------------------------------------------
  2411.  
  2412. [Function]
  2413. complex realpart &optional imagpart
  2414.  
  2415. The arguments must be non-complex numbers; a number is returned that has
  2416. realpart as its real part and imagpart as its imaginary part, possibly
  2417. converted according to the rule of floating-point contagion (thus both
  2418. components will be of the same type). If imagpart is not specified, then
  2419. (coerce 0 (type-of realpart)) is effectively used. Note that if both the
  2420. realpart and imagpart are rational and the imagpart is zero, then the result is
  2421. just the realpart because of the rule of canonical representation for complex
  2422. rationals. It follows that the result of complex is not always a complex
  2423. number; it may be simply a rational.
  2424.  
  2425. [Function]
  2426. realpart number
  2427. imagpart number
  2428.  
  2429. These return the real and imaginary parts of a complex number. If number is a
  2430. non-complex number, then realpart returns its argument number and imagpart
  2431. returns (* 0 number), which has the effect that the imaginary part of a
  2432. rational is 0 and that of a floating-point number is a floating-point zero of
  2433. the same format.
  2434.  
  2435. [change_begin]
  2436. A clever way to multiply a complex number z by i is to write
  2437.  
  2438. (complex (- (imagpart z)) (realpart z))
  2439.  
  2440. instead of (* z #c(0 1)). This cleverness is not always gratuitous; it may be
  2441. of particular importance in the presence of minus zero. For example, if we are
  2442. using IEEE standard floating-point arithmetic and z=4+0i, the result of the
  2443. clever expression is -0+4i, a true    rotation of z, whereas the result of (* z
  2444. #c(0 1)) is likely to be
  2445.  
  2446.  
  2447. (4+0i)(+0+i) = ((4)(+0)-(+0)(1))+((4)(1)+(+0)(+0)i
  2448.      = ((+0)-(+0))+((4)+(+0))i = +0+4i
  2449.  
  2450. which could land on the wrong side of a branch cut, for example.
  2451. [change_end]
  2452.  
  2453. -------------------------------------------------------------------------------
  2454.  
  2455. 12.7. Logical Operations on Numbers
  2456.  
  2457. The logical operations in this section require integers as arguments; it is an
  2458. error to supply a non-integer as an argument. The functions all treat integers
  2459. as if they were represented in two's-complement notation.
  2460.  
  2461. -------------------------------------------------------------------------------
  2462. Implementation note: Internally, of course, an implementation of Common Lisp
  2463. may or may not use a two's-complement representation. All that is necessary is
  2464. that the logical operations perform calculations so as to give this appearance
  2465. to the user.
  2466. -------------------------------------------------------------------------------
  2467.  
  2468. The logical operations provide a convenient way to represent an infinite vector
  2469. of bits. Let such a conceptual vector be indexed by the non-negative integers.
  2470. Then bit j is assigned a ``weight''   . Assume that only a finite number of
  2471. bits are 1's or only a finite number of bits are 0's. A vector with only a
  2472. finite number of one-bits is represented as the sum of the weights of the
  2473. one-bits, a positive integer. A vector with only a finite number of zero-bits
  2474. is represented as -1 minus the sum of the weights of the zero-bits, a negative
  2475. integer.
  2476.  
  2477. This method of using integers to represent bit-vectors can in turn be used to
  2478. represent sets. Suppose that some (possibly countably infinite) universe of
  2479. discourse for sets is mapped into the non-negative integers. Then a set can be
  2480. represented as a bit vector; an element is in the set if the bit whose index
  2481. corresponds to that element is a one-bit. In this way all finite sets can be
  2482. represented (by positive integers), as well as all sets whose complements are
  2483. finite (by negative integers). The functions logior, logand, and logxor defined
  2484. below then compute the union, intersection, and symmetric difference operations
  2485. on sets represented in this way.
  2486.  
  2487. [Function]
  2488. logior &rest integers
  2489.  
  2490. This returns the bit-wise logical inclusive or of its arguments. If no argument
  2491. is given, then the result is zero, which is an identity for this operation.
  2492.  
  2493. [Function]
  2494. logxor &rest integers
  2495.  
  2496. This returns the bit-wise logical exclusive or of its arguments. If no argument
  2497. is given, then the result is zero, which is an identity for this operation.
  2498.  
  2499. [Function]
  2500. logand &rest integers
  2501.  
  2502. This returns the bit-wise logical and of its arguments. If no argument is
  2503. given, then the result is -1, which is an identity for this operation.
  2504.  
  2505. [Function]
  2506. logeqv &rest integers
  2507.  
  2508. This returns the bit-wise logical equivalence (also known as exclusive nor) of
  2509. its arguments. If no argument is given, then the result is -1, which is an
  2510. identity for this operation.
  2511.  
  2512. [Function]
  2513. lognand integer1 integer2
  2514. lognor integer1 integer2
  2515. logandc1 integer1 integer2
  2516. logandc2 integer1 integer2
  2517. logorc1 integer1 integer2
  2518. logorc2 integer1 integer2
  2519.  
  2520. These are the other six non-trivial bit-wise logical operations on two
  2521. arguments. Because they are not associative, they take exactly two arguments
  2522. rather than any non-negative number of arguments.
  2523.  
  2524. (lognand n1 n2) == (lognot (logand n1 n2))
  2525. (lognor n1 n2) == (lognot (logior n1 n2))
  2526. (logandc1 n1 n2) == (logand (lognot n1) n2)
  2527. (logandc2 n1 n2) == (logand n1 (lognot n2))
  2528. (logorc1 n1 n2) == (logior (lognot n1) n2)
  2529. (logorc2 n1 n2) == (logior n1 (lognot n2))
  2530.  
  2531. The ten bit-wise logical operations on two integers are summarized in the
  2532. following table:
  2533.  
  2534. ----------------------------------------------------------------
  2535. integer1        0       0       1       1
  2536. integer2        0       1       0       1       Operation Name
  2537. ----------------------------------------------------------------
  2538. logand          0       0       0       1       and
  2539. logior          0       1       1       1       inclusive or
  2540. logxor          0       1       1       0       exclusive or
  2541. logeqv          1       0       0       1       equivalence (exclusive nor)
  2542. lognand         1       1       1       0       not-and
  2543. lognor          1       0       0       0       not-or
  2544. logandc1        0       1       0       0       and complement of integer1 with integer2
  2545. logandc2        0       0       1       0       and integer1 with complement of integer2
  2546. logorc1         1       1       0       1       or complement of integer1 with integer2
  2547. logorc2         1       0       1       1       or integer1 with complement of integer2
  2548.  
  2549. [Function]
  2550. boole op integer1 integer2
  2551.  
  2552. [Constant]
  2553. boole-clr 
  2554. boole-set 
  2555. boole-1 
  2556. boole-2 
  2557. boole-c1 
  2558. boole-c2 
  2559. boole-and 
  2560. boole-ior 
  2561. boole-xor 
  2562. boole-eqv 
  2563. boole-nand 
  2564. boole-nor 
  2565. boole-andc1 
  2566. boole-andc2 
  2567. boole-orc1 
  2568. boole-orc2 
  2569.  
  2570. The function boole takes an operation op and two integers, and returns an
  2571. integer produced by performing the logical operation specified by op on the two
  2572. integers. The precise values of the sixteen constants are
  2573. implementation-dependent, but they are suitable for use as the first argument
  2574. to boole:
  2575.  
  2576. ----------------------------------------------------------------
  2577. integer1        0       0       1       1
  2578. integer2        0       1       0       1       Operation Performed
  2579. ----------------------------------------------------------------
  2580. boole-clr       0       0       0       0       always 0
  2581. boole-set       1       1       1       1       always 1
  2582. boole-1         0       0       1       1       integer1
  2583. boole-2         0       1       0       1       integer2
  2584. boole-c1        1       1       0       0       complement of integer1
  2585. boole-c2        1       0       1       0       complement of integer2
  2586. boole-and       0       0       0       1       and
  2587. boole-ior       0       1       1       1       inclusive or
  2588. boole-xor       0       1       1       0       exclusive or
  2589. boole-eqv       1       0       0       1       equivalence (exclusive nor)
  2590. boole-nand      1       1       1       0       not-and
  2591. boole-nor       1       0       0       0       not-or
  2592. boole-andc1     0       1       0       0       and complement of integer1 with integer2
  2593. boole-andc2     0       0       1       0       and integer1 with complement of integer2
  2594. boole-orc1      1       1       0       1       or complement of integer1 with integer2
  2595. boole-orc2      1       0       1       1       or integer1 with complement of integer2
  2596.  
  2597. boole can therefore compute all sixteen logical functions on two arguments. In
  2598. general,
  2599.  
  2600. (boole boole-and x y) == (logand x y)
  2601.  
  2602. and the latter is more perspicuous. However, boole is useful when it is
  2603. necessary to parameterize a procedure so that it can use one of several logical
  2604. operations.
  2605.  
  2606. [Function]
  2607. lognot integer
  2608.  
  2609. This returns the bit-wise logical not of its argument. Every bit of the result
  2610. is the complement of the corresponding bit in the argument.
  2611.  
  2612. (logbitp j (lognot x)) == (not (logbitp j x))
  2613.  
  2614. [Function]
  2615. logtest integer1 integer2
  2616.  
  2617. logtest is a predicate that is true if any of the bits designated by the 1's in
  2618. integer1 are 1's in integer2.
  2619.  
  2620. (logtest x y) == (not (zerop (logand x y)))
  2621.  
  2622. [Function]
  2623. logbitp index integer
  2624.  
  2625. logbitp is true if the bit in integer whose index is index (that is, its weight
  2626. is   ) is a one-bit; otherwise it is false. For example:
  2627.  
  2628. (logbitp 2 6) is true
  2629. (logbitp 0 6) is false
  2630. (logbitp k n) == (ldb-test (byte 1 k) n)
  2631.  
  2632. [change_begin]
  2633. X3J13 voted in January 1989 (ARGUMENTS-UNDERSPECIFIED)   to clarify that the
  2634. index must be a non-negative integer.
  2635. [change_end]
  2636.  
  2637. [Function]
  2638. ash integer count
  2639.  
  2640. This function shifts integer arithmetically left by count bit positions if
  2641. count is positive, or right by -count bit positions if count is negative. The
  2642. sign of the result is always the same as the sign of integer.
  2643.  
  2644. Mathematically speaking, this operation performs the computation   ).
  2645.  
  2646. Logically, this moves all of the bits in integer to the left, adding zero-bits
  2647. at the bottom, or moves them to the right, discarding bits. (In this context
  2648. the question of what gets shifted in on the left is irrelevant; integers,
  2649. viewed as strings of bits, are ``half-infinite,'' that is, conceptually extend
  2650. infinitely far to the left.) For example:
  2651.  
  2652. (logbitp j (ash n k)) == (and (>= j k) (logbitp (- j k) n))
  2653.  
  2654. [Function]
  2655. logcount integer
  2656.  
  2657. The number of bits in integer is determined and returned. If integer is
  2658. positive, the 1-bits in its binary representation are counted. If integer is
  2659. negative, the 0-bits in its two's-complement binary representation are counted.
  2660. The result is always a non-negative integer. For example:
  2661.  
  2662. (logcount 13) => 3      ;Binary representation is ...0001101
  2663. (logcount -13) => 2     ;Binary representation is ...1110011
  2664. (logcount 30) => 4      ;Binary representation is ...0011110
  2665. (logcount -30) => 4     ;Binary representation is ...1100010
  2666.  
  2667. The following identity always holds:
  2668.  
  2669. (logcount x) == (logcount (- (+ x 1)))
  2670.              == (logcount (lognot x))
  2671.  
  2672. [Function]
  2673. integer-length integer
  2674.  
  2675. This function performs the computation
  2676.  
  2677.  
  2678.  
  2679. This is useful in two different ways. First, if integer is non-negative, then
  2680. its value can be represented in unsigned binary form in a field whose width in
  2681. bits is no smaller than (integer-length integer). Second, regardless of the
  2682. sign of integer, its value can be represented in signed binary two's-complement
  2683. form in a field whose width in bits is no smaller than (+ (integer-length
  2684. integer) 1). For example:
  2685.  
  2686. (integer-length 0) => 0
  2687. (integer-length 1) => 1
  2688. (integer-length 3) => 2
  2689. (integer-length 4) => 3
  2690. (integer-length 7) => 3
  2691. (integer-length -1) => 0
  2692. (integer-length -4) => 2
  2693. (integer-length -7) => 3
  2694. (integer-length -8) => 3
  2695.  
  2696. -------------------------------------------------------------------------------
  2697. Compatibility note: This function is similar to the MacLisp function haulong.
  2698. One may define haulong as
  2699.  
  2700. (haulong x) == (integer-length (abs x))
  2701.  
  2702. -------------------------------------------------------------------------------
  2703.  
  2704. -------------------------------------------------------------------------------
  2705.  
  2706. 12.8. Byte Manipulation Functions
  2707.  
  2708. Several functions are provided for dealing with an arbitrary-width field of
  2709. contiguous bits appearing anywhere in an integer. Such a contiguous set of bits
  2710. is called a byte. Here the term byte does not imply some fixed number of bits
  2711. (such as eight), rather a field of arbitrary and user-specifiable width.
  2712.  
  2713. The byte-manipulation functions use objects called byte specifiers to designate
  2714. a specific byte position within an integer. The representation of a byte
  2715. specifier is implementation-dependent; in particular, it may or may not be a
  2716. number. It is sufficient to know that the function byte will construct one, and
  2717. that the byte-manipulation functions will accept them. The function byte
  2718. accepts two integers representing the position and size of the byte and returns
  2719. a byte specifier. Such a specifier designates a byte whose width is size and
  2720. whose bits have weights    through   .
  2721.  
  2722. [Function]
  2723. byte size position
  2724.  
  2725. byte takes two integers representing the size and position of a byte and
  2726. returns a byte specifier suitable for use as an argument to byte-manipulation
  2727. functions.
  2728.  
  2729. [Function]
  2730. byte-size bytespec
  2731. byte-position bytespec
  2732.  
  2733. Given a byte specifier, byte-size returns the size specified as an integer;
  2734. byte-position similarly returns the position. For example:
  2735.  
  2736. (byte-size (byte j k)) == j
  2737. (byte-position (byte j k)) == k
  2738.  
  2739. [Function]
  2740. ldb bytespec integer
  2741.  
  2742. bytespec specifies a byte of integer to be extracted. The result is returned as
  2743. a non-negative integer. For example:
  2744.  
  2745. (logbitp j (ldb (byte s p) n)) == (and (< j s) (logbitp (+ j p) n))
  2746.  
  2747. The name of the function ldb means ``load byte.''
  2748.  
  2749. -------------------------------------------------------------------------------
  2750. Compatibility note: The MacLisp function haipart can be implemented in terms of
  2751. ldb as follows:
  2752.  
  2753. (defun haipart (integer count)
  2754.   (let ((x (abs integer)))
  2755.     (if (minusp count)
  2756.         (ldb (byte (- count) 0) x)
  2757.         (ldb (byte count (max 0 (- (integer-length x) count)))
  2758.              x))))
  2759.  
  2760. -------------------------------------------------------------------------------
  2761.  
  2762. If the argument integer is specified by a form that is a place form acceptable
  2763. to setf, then setf may be used with ldb to modify a byte within the integer
  2764. that is stored in that place. The effect is to perform a dpb operation and then
  2765. store the result back into the place.
  2766.  
  2767. [Function]
  2768. ldb-test bytespec integer
  2769.  
  2770. ldb-test is a predicate that is true if any of the bits designated by the byte
  2771. specifier bytespec are 1's in integer; that is, it is true if the designated
  2772. field is non-zero.
  2773.  
  2774. (ldb-test bytespec n) == (not (zerop (ldb bytespec n)))
  2775.  
  2776. [Function]
  2777. mask-field bytespec integer
  2778.  
  2779. This is similar to ldb; however, the result contains the specified byte of
  2780. integer in the position specified by bytespec, rather than in position 0 as
  2781. with ldb. The result therefore agrees with integer in the byte specified but
  2782. has zero-bits everywhere else. For example:
  2783.  
  2784. (ldb bs (mask-field bs n)) == (ldb bs n)
  2785.  
  2786. (logbitp j (mask-field (byte s p) n))
  2787.    == (and (>= j p) (< j (+ p s)) (logbitp j n))
  2788.  
  2789. (mask-field bs n) == (logand n (dpb -1 bs 0))
  2790.  
  2791. If the argument integer is specified by a form that is a place form acceptable
  2792. to setf, then setf may be used with mask-field to modify a byte within the
  2793. integer that is stored in that place. The effect is to perform a deposit-field
  2794. operation and then store the result back into the place.
  2795.  
  2796. [Function]
  2797. dpb newbyte bytespec integer
  2798.  
  2799. This returns a number that is the same as integer except in the bits specified
  2800. by bytespec. Let s be the size specified by bytespec; then the low s bits of
  2801. newbyte appear in the result in the byte specified by bytespec. The integer
  2802. newbyte is therefore interpreted as being right-justified, as if it were the
  2803. result of ldb. For example:
  2804.  
  2805. (logbitp j (dpb m (byte s p) n))
  2806.   == (if (and (>= j p) (< j (+ p s)))
  2807.          (logbitp (- j p) m)
  2808.          (logbitp j n))
  2809.  
  2810. The name of the function dpb means ``deposit byte.''
  2811.  
  2812. [Function]
  2813. deposit-field newbyte bytespec integer
  2814.  
  2815. This function is to mask-field as dpb is to ldb. The result is an integer that
  2816. contains the bits of newbyte within the byte specified by bytespec, and
  2817. elsewhere contains the bits of integer. For example:
  2818.  
  2819. (logbitp j (deposit-field m (byte s p) n))
  2820.    == (if (and (>= j p) (< j (+ p s)))
  2821.           (logbitp j m)
  2822.           (logbitp j n))
  2823.  
  2824. -------------------------------------------------------------------------------
  2825. Implementation note: If the bytespec is a constant, one may of course
  2826. construct, at compile time, an equivalent mask m, for example by computing
  2827. (deposit-field -1 bytespec 0). Given this mask m, one may then compute
  2828.  
  2829. (deposit-field newbyte bytespec integer)
  2830.  
  2831. by computing
  2832.  
  2833. (logior (logand newbyte m) (logand integer (lognot m)))
  2834.  
  2835. where the result of (lognot m) can of course also be computed at compile time.
  2836. However, the following expression may also be used and may require fewer
  2837. temporary registers in some situations:
  2838.  
  2839. (logxor integer (logand m (logxor integer newbyte)))
  2840.  
  2841. A related, though possibly less useful, trick is that
  2842.  
  2843. (let ((z (logand (logxor x y) m)))
  2844.   (setq x (logxor z x))
  2845.   (setq y (logxor z y)))
  2846.  
  2847. interchanges those bits of x and y for which the mask m is 1, and leaves alone
  2848. those bits of x and y for which m is 0.
  2849. -------------------------------------------------------------------------------
  2850.  
  2851. -------------------------------------------------------------------------------
  2852.  
  2853. 12.9. Random Numbers
  2854.  
  2855. The Common Lisp facility for generating pseudo-random numbers has been
  2856. carefully defined to make its use reasonably portable. While two
  2857. implementations may produce different series of pseudo-random numbers, the
  2858. distribution of values should be relatively independent of such
  2859. machine-dependent aspects as word size.
  2860.  
  2861. [Function]
  2862. random number &optional state
  2863.  
  2864. (random n) accepts a positive number n and returns a number of the same kind
  2865. between zero (inclusive) and n (exclusive). The number n may be an integer or a
  2866. floating-point number. An approximately uniform choice distribution is used. If
  2867. n is an integer, each of the possible results occurs with (approximate)
  2868. probability
  2869. 1/n. (The qualifier ``approximate'' is used because of implementation
  2870. considerations; in practice, the deviation from uniformity should be quite
  2871. small.)
  2872.  
  2873. The argument state must be an object of type random-state; it defaults to the
  2874. value of the variable *random-state*. This object is used to maintain the state
  2875. of the pseudo-random-number generator and is altered as a side effect of the
  2876. random operation.
  2877.  
  2878. -------------------------------------------------------------------------------
  2879. Compatibility note: random of zero arguments as defined in MacLisp has been
  2880. omitted because its value is too implementation-dependent (limited by fixnum
  2881. range).
  2882. -------------------------------------------------------------------------------
  2883. Implementation note: In general, even if random of zero arguments were defined
  2884. as in MacLisp, it is not adequate to define (random n) for integral n to be
  2885. simply (mod (random) n); this fails to be uniformly distributed if n is larger
  2886. than the largest number produced by random, or even if n merely approaches this
  2887. number. This is another reason for omitting random of zero arguments in Common
  2888. Lisp. Assuming that the underlying mechanism produces ``random bits'' (possibly
  2889. in chunks such as fixnums), the best approach is to produce enough random bits
  2890. to construct an integer k some number d of bits larger than (integer-length n)
  2891. (see integer-length), and then compute (mod k n). The quantity d should be at
  2892. least 7, and preferably 10 or more.
  2893.  
  2894. To produce random floating-point numbers in the half-open range [A, B),
  2895. accepted practice (as determined by a look through the Collected Algorithms
  2896. from the ACM, particularly algorithms 133, 266, 294, and 370) is to compute
  2897. X * (B - A) + A, where X is a floating-point number uniformly distributed over
  2898. [0.0, 1.0) and computed by calculating a random integer N in the range [0, M)
  2899. (typically by a multiplicative-congruential or linear-congruential method mod
  2900. M) and then setting X = N/M. See also [27]. If one takes   , where f is the
  2901. length of the significand of a floating-point number (and it is in fact common
  2902. to choose M to be a power of 2), then this method is equivalent to the
  2903. following assembly-language-level procedure. Assume the representation has no
  2904. hidden bit. Take a floating-point 0.5, and clobber its entire significand with
  2905. random bits. Normalize the result if necessary.
  2906.  
  2907. For example, on the DEC PDP-10, assume that accumulator T is completely random
  2908. (all 36 bits are random). Then the code sequence
  2909.  
  2910. LSH T,-9                 ;Clear high 9 bits; low 27 are random
  2911. FSC T,128.               ;Install exponent and normalize
  2912.  
  2913. will produce in T a random floating-point number uniformly distributed over
  2914. [0.0, 1.0). (Instead of the LSH instruction, one could do
  2915.  
  2916. TLZ T,777000             ;That's 777000 octal
  2917.  
  2918. but if the 36 random bits came from a congruential random-number generator, the
  2919. high-order bits tend to be ``more random'' than the low-order ones, and so the
  2920. LSH would be better for uniform distribution. Ideally all the bits would be the
  2921. result of high-quality randomness.)
  2922.  
  2923. With a hidden-bit representation, normalization is not a problem, but dealing
  2924. with the hidden bit is. The method can be adapted as follows. Take a
  2925. floating-point 1.0 and clobber the explicit significand bits with random bits;
  2926. this produces a random floating-point number in the range [1.0, 2.0). Then
  2927. simply subtract 1.0. In effect, we let the hidden bit creep in and then
  2928. subtract it away again.
  2929.  
  2930. For example, on the DEC VAX, assume that register T is completely random (but a
  2931. little less random than on the PDP-10, as it has only 32 random bits). Then the
  2932. code sequence
  2933.  
  2934. INSV #^X81,#7,#9,T     ;Install correct sign bit and exponent
  2935. SUBF #^F1.0,T          ;Subtract 1.0
  2936.  
  2937. will produce in T a random floating-point number uniformly distributed over
  2938. [0.0, 1.0). Again, if the low-order bits are not random enough, then the
  2939. instruction
  2940.  
  2941. ROTL #7,T
  2942.  
  2943. should be performed first.
  2944.  
  2945. Implementors may wish to consult reference [41] for a discussion of some
  2946. efficient methods of generating pseudo-random numbers.
  2947. -------------------------------------------------------------------------------
  2948.  
  2949. [Variable]
  2950. *random-state*
  2951.  
  2952. This variable holds a data structure, an object of type random-state, that
  2953. encodes the internal state of the random-number generator that random uses by
  2954. default. The nature of this data structure is implementation-dependent. It may
  2955. be printed out and successfully read back in, but may or may not function
  2956. correctly as a random-number state object in another implementation. A call to
  2957. random will perform a side effect on this data structure. Lambda-binding this
  2958. variable to a different random-number state object will correctly save and
  2959. restore the old state object.
  2960.  
  2961. [Function]
  2962. make-random-state &optional state
  2963.  
  2964. This function returns a new object of type random-state, suitable for use as
  2965. the value of the variable *random-state*. If state is nil or omitted,
  2966. make-random-state returns a copy of the current random-number state object (the
  2967. value of the variable *random-state*). If state is a state object, a copy of
  2968. that state object is returned. If state is t, then a new state object is
  2969. returned that has been ``randomly'' initialized by some means (such as by a
  2970. time-of-day clock).
  2971.  
  2972. -------------------------------------------------------------------------------
  2973. Rationale: Common Lisp purposely provides no way to initialize a random-state
  2974. object from a user-specified ``seed.'' The reason for this is that the number
  2975. of bits of state information in a random-state object may vary widely from one
  2976. implementation to another, and there is no simple way to guarantee that any
  2977. user-specified seed value will be ``random enough.'' Instead, the
  2978. initialization of random-state objects is left to the implementor in the case
  2979. where the argument t is given to make-random-state.
  2980.  
  2981. To handle the common situation of executing the same program many times in a
  2982. reproducible manner, where that program uses random, the following procedure
  2983. may be used:
  2984.  
  2985.   1.  Evaluate (make-random-state t) to create a random-state object.
  2986.  
  2987.   2.  Write that object to a file, using print, for later use.
  2988.  
  2989.   3.  Whenever the program is to be run, first use read to create a copy of the
  2990.      random-state object from the printed representation in the file. Then use
  2991.      the random-state object newly created by the read operation to initialize
  2992.      the random-number generator for the program.
  2993.  
  2994. It is for the sake of this procedure for reproducible execution that
  2995. implementations are required to provide a read/print syntax for objects of type
  2996. random-state.
  2997.  
  2998. It is also possible to make copies of a random-state object directly without
  2999. going through the print/read process, simply by using the make-random-state
  3000. function to copy the object; this allows the same sequence of random numbers to
  3001. be generated many times within a single program.
  3002. -------------------------------------------------------------------------------
  3003. Implementation note: A recommended way to implement the type random-state is
  3004. effectively to use the machinery for defstruct. The usual structure syntax may
  3005. then be used for printing random-state objects; one might look something like
  3006.  
  3007. #S(RANDOM-STATE DATA #(14 49 98436589 786345 8734658324 ...))
  3008.  
  3009. where the components are of course completely implementation-dependent.
  3010. -------------------------------------------------------------------------------
  3011.  
  3012. [Function]
  3013. random-state-p object
  3014.  
  3015. random-state-p is true if its argument is a random-state object, and otherwise
  3016. is false.
  3017.  
  3018. (random-state-p x) == (typep x 'random-state)
  3019.  
  3020. -------------------------------------------------------------------------------
  3021.  
  3022. 12.10 Implementation Parameters
  3023.  
  3024. The values of the named constants defined in this section are
  3025. implementation-dependent. They may be useful for parameterizing code in some
  3026. situations.
  3027.  
  3028. [Constant]
  3029. most-positive-fixnum 
  3030. most-negative-fixnum 
  3031.  
  3032. The value of most-positive-fixnum is that fixnum closest in value to positive
  3033. infinity provided by the implementation.
  3034.  
  3035. The value of most-negative-fixnum is that fixnum closest in value to negative
  3036. infinity provided by the implementation.
  3037.  
  3038. [change_begin]
  3039. X3J13 voted in January 1989 (FIXNUM-NON-PORTABLE)   to specify that fixnum must
  3040. be a supertype of the type (signed-byte 16), and additionally that the value of
  3041. array-dimension-limit must be a fixnum. This implies that the value of
  3042. most-negative-fixnum must be less than or equal to   , and the value of
  3043. most-positive-fixnum must be greater than or equal to both    and the value of
  3044. array-dimension-limit.
  3045. [change_end]
  3046.  
  3047. [Constant]
  3048. most-positive-short-float 
  3049. least-positive-short-float 
  3050. least-negative-short-float 
  3051. most-negative-short-float 
  3052.  
  3053. The value of most-positive-short-float is that short-format floating-point
  3054. number closest in value to (but not equal to) positive infinity provided by the
  3055. implementation.
  3056.  
  3057. The value of least-positive-short-float is that positive short-format
  3058. floating-point number closest in value to (but not equal to) zero provided by
  3059. the implementation.
  3060.  
  3061. The value of least-negative-short-float is that negative short-format
  3062. floating-point number closest in value to (but not equal to) zero provided by
  3063. the implementation. (Note that even if an implementation supports minus zero as
  3064. a distinct short floating-point value, least-negative-short-float must not be
  3065. minus zero.)
  3066.  
  3067. [change_begin]
  3068. X3J13 voted in June 1989 (FLOAT-UNDERFLOW)   to clarify that these definitions
  3069. are to be taken quite literally. In implementations that support denormalized
  3070. numbers, the values of least-positive-short-float and
  3071. least-negative-short-float may be denormalized.
  3072. [change_end]
  3073.  
  3074. The value of most-negative-short-float is that short-format floating-point
  3075. number closest in value to (but not equal to) negative infinity provided by the
  3076. implementation.
  3077.  
  3078. [Constant]
  3079. most-positive-single-float 
  3080. least-positive-single-float 
  3081. least-negative-single-float 
  3082. most-negative-single-float 
  3083. most-positive-double-float 
  3084. least-positive-double-float 
  3085. least-negative-double-float 
  3086. most-negative-double-float 
  3087. most-positive-long-float 
  3088. least-positive-long-float 
  3089. least-negative-long-float 
  3090. most-negative-long-float 
  3091.  
  3092. These are analogous to the constants defined above for short-format
  3093. floating-point numbers.
  3094.  
  3095. [change_begin]
  3096. [Constant]
  3097. least-positive-normalized-short-float 
  3098. least-negative-normalized-short-float 
  3099.  
  3100. X3J13 voted in June 1989 (FLOAT-UNDERFLOW)   to add these constants to the
  3101. language.
  3102.  
  3103. The value of least-positive-normalized-short-float is that positive normalized
  3104. short-format floating-point number closest in value to (but not equal to) zero
  3105. provided by the implementation. In implementations that do not support
  3106. denormalized numbers this may be the same as the value of
  3107. least-positive-short-float.
  3108.  
  3109. The value of least-negative-normalized-short-float is that negative normalized
  3110. short-format floating-point number closest in value to (but not equal to) zero
  3111. provided by the implementation. (Note that even if an implementation supports
  3112. minus zero as a distinct short floating-point value,
  3113. least-negative-normalized-short-float must not be minus zero.) In
  3114. implementations that do not support denormalized numbers this may be the same
  3115. as the value of least-positive-short-float.
  3116.  
  3117. [Constant]
  3118. least-positive-normalized-single-float 
  3119. least-negative-normalized-single-float 
  3120. least-positive-normalized-double-float 
  3121. least-negative-normalized-double-float 
  3122. least-positive-normalized-long-float 
  3123. least-negative-normalized-long-float 
  3124.  
  3125. These are analogous to the constants defined above for short-format
  3126. floating-point numbers.
  3127. [change_end]
  3128.  
  3129. [Constant]
  3130. short-float-epsilon 
  3131. single-float-epsilon 
  3132. double-float-epsilon 
  3133. long-float-epsilon 
  3134.  
  3135. These constants have as value, for each floating-point format, the smallest
  3136. positive floating-point number e of that format such that the expression
  3137.  
  3138. (not (= (float 1 e) (+ (float 1 e) e)))
  3139.  
  3140. is true when actually evaluated.
  3141.  
  3142. [Constant]
  3143. short-float-negative-epsilon 
  3144. single-float-negative-epsilon 
  3145. double-float-negative-epsilon 
  3146. long-float-negative-epsilon 
  3147.  
  3148. These constants have as value, for each floating-point format, the smallest
  3149. positive floating-point number e of that format such that the expression
  3150.  
  3151. (not (= (float 1 e) (- (float 1 e) e)))
  3152.  
  3153. is true when actually evaluated.
  3154.  
  3155. -------------------------------------------------------------------------------
  3156.  
  3157.  
  3158.  
  3159.  
  3160.  
  3161.  
  3162.